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Preface 


This  report  contains  reprints  of  four  papers  related  to  use  of  the  generalized  self-consistent 
method  (GSCM)  for  determining  the  homogenized  constitutive  response  of  microcracked  media 
for  the  development  of  multiscale  constitutive  models.  In  “The  Effect  of  Crack  Face  Contact  on 
the  Anisotropic  Effective  Moduli  of  Microcrack  Damaged  Media,”1  the  GSCM  is  used  in 
conjunction  with  a  finite-element  method  to  determine  the  anisotropic  effective  moduli  of  a 
medium  containing  damage  consisting  of  microcracks  with  an  arbitrary  degree  of  alignment. 

The  moduli  of  the  medium  subjected  to  tension,  compression,  and  an  initially  stress-free  state  are 
evaluated  and  shown  to  be  significantly  different,  affecting  the  wave  speed  (illustrated  using 
slowness  surfaces)  in  the  damaged  medium.  In  “An  Effective  Medium  Model  for  Elastic  Waves 
in  Microcrack  Damaged  Media,”2  direct  numerical  simulations  of  waves  traveling  in  microcrack 
damaged  media  are  compared  to  results  using  a  homogenized  effective  medium  calculation.  For 
incident  waves  with  large  wavelengths  (1/ka),  where  k  is  the  wavenumber  and  a  is  the  half-crack 
length,  the  scattered  elastic  energy  approaches  0  and  the  wave  does  not  “see”  the  obstacle;  for 
crack  systems  modeled  using  the  finite-element  approach,  this  occurs  for  1/ka  >  60  for  media  in 
tension  and  1/ka  >  10  for  those  in  compression.  In  “Anisotropic  Effective  Moduli  of 
Microcracked  Materials  Under  Antiplane  Loading,”3  the  anisotropic  effective  moduli  of  a 
cracked  solid  subjected  to  antiplane  shear  deformation  are  analytically  determined.  When  the 
undamaged  solid  is  isotropic,  the  GSCM  can  be  realized  exactly;  however,  when  the  undamaged 
solid  is  anisotropic,  coupled  nonlinear  equations  for  the  unknown  effective  moduli  are 
determined  through  numerical  iteration.  Finally,  in  “On  the  Effective  Electroelastic  Properties  of 
Microcracked  Generally  Anisotropic  Solids,”4  concise  expressions  are  derived  for  the  effective 
electroelastic  properties  of  a  piezoelectric  solid  containing  insulating,  permeable,  or  conducting 
microcracks. 
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Abstract 

The  generalized  self-consistent  method  (GSCM)  in  conjunction  with  a  computational  finite  element  method  is  used  to 
calculate  the  anisotropic  effective  moduli  of  a  medium  containing  damage  consisting  of  microcracks  with  an  arbitrary 
degree  of  alignment.  Since  cracks  respond  differently  under  different  external  loads,  the  moduli  of  the  medium  subjected 
to  tension,  compression  and  an  initially  stress-free  state  are  evaluated  and  shown  to  be  significantly  different,  which  will 
further  affect  the  wave  speed  inside  the  damaged  media.  There  are  four  independent  material  moduli  for  a  2-D  plane  stress 
orthotropic  medium  in  tension  or  compression,  and  seven  independent  material  moduli  for  a  2-D  plane  stress  orthotropic 
cracked  medium,  which  is  initially  stress  free.  When  friction  exists,  it  further  changes  the  effective  moduli.  Numerical  meth¬ 
ods  are  used  to  take  into  account  crack  face  contact  and  friction.  The  wave  slowness  profiles  for  microcrack  damaged 
media  are  plotted  using  the  predicted  effective  material  moduli. 

©  2006  Elsevier  Ltd.  All  rights  reserved. 

Keywords:  Anisotropic  damage;  Effective  moduli;  Generalized  self-consistent  method;  Microcracks 


1.  Introduction 


Microcracking  is  a  major  cause  of  damage  in  brittle  material  systems,  which  reduces  the  material’s  effective 
moduli  by  allowing  an  increase  in  the  large-scale,  average  deformation  resulting  from  a  given  external  load.  As 
the  number  or  size  of  the  microcracks  increase,  there  is  a  commensurate  reduction  in  the  effective  stiffness, 
until  the  damage  reaches  a  critical  value  causing  the  material  to  fail.  A  convenient  way  to  quantify  the  crack 
density  in  2-D  is  the  parameter  tj  [1]. 


(i) 
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where  Mis  the  number  of  cracks  per  unit  area  A,  and  c  is  the  average  half  crack  length.  In  the  general  case,  an 
ensemble  of  microcracks  can  be  characterized  by  a  distribution  function,  which  includes  the  average  orienta¬ 
tion  of  the  cracks  and  a  description  of  how  the  crack  angles  are  distributed  around  the  average.  This  will  result 
in  anisotropic  effective  moduli,  even  when  the  undamaged  material  behaves  isotropically.  Since  microcracks 
respond  differently  under  tension  and  compression,  damage  growth  and  the  resulting  anisotropic  effective 
moduli  of  microcracked  media  will  differ  depending  upon  the  magnitude  and  direction  of  the  boundary  loads. 

The  crack-induced  changes  in  the  anisotropic  effective  moduli  will  have  a  profound  effect  on  the  overall 
structural  response  to  loads.  The  prediction  of  these  changes  will  help  in  the  understanding  of  the  anisotropic 
damage  evolution  process  in  brittle  materials  such  as  concrete,  rock  and  ceramics.  Indeed,  continuum  damage 
finite  element  modeling  of  the  dynamic  response  of  such  brittle  materials  often  requires,  at  the  very  least,  an 
orthotropic  description  of  damage  in  terms  of  three  orthogonal  arrays  of  interacting  cracks  [2].  The  aniso¬ 
tropic  response,  is  in  part,  a  load-induced  effect  [3],  but  can  also  be  attributed  to  a  pre-existing  preferred  ori¬ 
entation  of  microcracks.  In  addition,  the  mechanical  properties  of  microcracked  media  have  a  significant  effect 
on  wave  propagation  even  in  cases  where  the  wavelengths  are  large  enough  to  ignore  wave  reflection  from  the 
crack  faces  (see  e.g.  [4-6]). 

Sayers  et  al.  [7,8]  used  the  effective  moduli  of  the  damaged  media  with  perfectly  aligned  or  randomly  ori¬ 
ented  cracks  to  study  the  elastic  wave  propagation.  Schubnel  and  Gueguen  [9]  and  Gueguen  and  Schubnel  [10] 
calculated  the  elastic  wave  velocities  and  permeabilities  in  microcracked  rocks.  However,  because  the  interac¬ 
tions  of  arbitrary  orientations  of  cracks  were  not  fully  characterized,  these  authors  only  focused  on  the  prob¬ 
lems  of  aligned,  perfectly  random  or  specially  distributed  cracks,  and  did  not  consider  the  difference  between 
the  tensile  and  compressive  behavior  of  microcracks.  Therefore,  further  work  is  needed  to  understand  the 
complex  relationships  among  microcrack  orientation  distributions,  external  loads  and  the  resulting  effective 
moduli. 

Many  researchers  have  proposed  methods  to  calculate  the  effective  moduli  for  solids  with  inhomogeneities, 
voids  or  cracks.  Without  attempting  to  review  the  vast  literature  on  the  subject,  we  will  focus  on  the  studies 
that  had  the  largest  direct  impact  on  the  current  work;  the  interested  reader  is  referred  to  a  comprehensive 
review  of  the  literature  provided  by  Kachanov  [11,12].  The  first  solution  that  considered  the  elastic  interaction 
of  elliptically  shaped  cracks  was  by  Budiansky  and  O’Connell  [13].  In  that  landmark  paper,  they  used  the  self- 
consistent  method  (SCM)  to  calculate  the  3-D  elastic  moduli  of  a  cracked  solid,  where  the  distribution  of 
microcracks  was  assumed  to  be  isotropic  and  homogeneous.  In  the  SCM,  a  crack  is  embedded  directly  into 
an  effective  medium  and  the  energy  associated  with  opening  the  crack  under  a  given  load  is  calculated.  By 
equating  the  strain  energy  associated  with  the  opening  of  a  distribution  of  such  cracks,  with  the  difference 
in  strain  energy  between  the  undamaged  and  the  effective  medium,  the  effective  moduli  for  the  cracked  med¬ 
ium  can  be  calculated. 

Gottesman  et  al.  [14]  used  variational  principles  to  calculate  the  effective  elastic  moduli  for  materials  with 
both  interacting  and  non-interacting  parallel  cracks  using  the  SCM.  These  two  solutions  provided  upper 
and  lower  bounds  for  the  effective  moduli  of  a  medium  with  aligned  cracks,  resulting  in  an  orthotropic  effective 
medium  as  opposed  to  the  randomly  oriented  cracks,  which  results  in  an  isotropic  effective  medium.  Horii  and 
Nemat-Nasser  [3]  extended  the  Budiansky  and  O’Connell  [13]  solution,  to  account  for  crack  closure  and  crack 
face  friction  effects.  They  considered  load  cases  in  which  some  cracks  opened  and  some  closed,  and  showed  that 
when  crack  closure  and  friction  are  considered,  the  overall  material  response  depends  on  load  history. 

Kachanov  [11]  provided  a  comprehensive  review  of  earlier  results,  and  compared  them  to  analytical  solu¬ 
tions  he  derived  for  non-interacting  cracks.  He  noted  that  there  are  several  drawbacks  in  the  use  of  the  SCM 
for  determining  effective  moduli  of  cracked  solids:  (1)  Because  this  method  inserts  the  cracks  directly  into  the 
effective  media,  the  SCM  over  estimates  the  interaction  between  the  cracks.  (2)  The  SCM  has  a  limiting  value 
for  crack  density  for  randomly  distributed  cracks:  it  gives  zero  stiffness  when  the  2-D  crack  density  r\  is  larger 
than  l/7i. 

In  order  to  overcome  these  difficulties,  several  researchers  have  developed  different  methods  for  determining 
the  effective  moduli  of  microcracked  media.  Benveniste  [15]  extended  the  effective  field  or  Mori-Tanaka 
method  (MTM)  [16]  to  calculate  the  effective  moduli  of  2-D  cracked  media.  This  MTM  has  some  similarities 
to  the  SCM,  but  the  energy  is  calculated  based  on  a  model  of  a  crack  in  an  undamaged  medium  subjected  to 
an  effective  far-field  stress  or  strain.  The  MTM  predicts  higher  moduli  than  the  SCM  and  does  not  predict  a 
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limiting  value  for  the  crack  density.  Hashin  [17]  developed  a  differential  scheme  method  (DSM)  for  elastic 
properties  of  a  cracked  material  by  following  a  limiting  process  on  the  DSM  for  porous  materials.  Kachanov 
[11]  points  out  that  the  DSM  is  an  incremental  version  of  the  SCM,  which  does  not  have  a  limiting  crack  den¬ 
sity  value,  but  still  under-predicts  the  effective  moduli,  although  not  as  severely  as  the  SCM.  At  the  same  time, 
he  notes  that  the  DSM  is  path  dependent,  since  it  depends  on  the  order  in  which  the  cracks  are  added  to  the 
formulation. 

Aboudi  and  Benveniste  [1]  used  a  generalized  self-consistent  method  (GSCM)  to  evaluate  the  2-D  effective 
moduli  of  a  material  with  randomly  oriented  cracks.  This  method  assumes  that  each  crack  resides  in  a  region 
of  undamaged  material,  which  in  turn  is  surrounded  by  the  effective  medium.  The  general  solution  for  a  crack 
inside  a  circular  inclusion  was  used  as  the  basis  for  this  work  (see  Fig.  1).  By  using  this  assumption,  the  GSCM 
gives  higher  effective  moduli  than  the  SCM  and  gives  non-zero  stiffness  at  crack  density  r\  =  \/n.  However, 
because  of  the  use  of  a  circular  inclusion,  the  crack  tip  touches  the  inclusion-effective  medium  interface  at 
rj  =  l/7i.  Therefore,  this  method  can  only  be  used  for  crack  densities  less  than  l/n  and  becomes  increasingly 
less  accurate  as  the  crack  tip  approaches  the  interface  r\  — ►  l/n.  Huang  et  al.  [18]  provided  a  crack-elliptical 
matrix-composite  model  for  the  GSCM;  in  this  improved  GSCM,  an  elliptical  inclusion  is  used  instead  of  a 
circular  inclusion,  which  extends  the  crack  density  limit  seen  in  Aboudi  and  Benveniste  [1]  to  a  higher  limiting 
value.  The  solution,  however,  is  still  only  presented  for  a  randomly  oriented  distribution  of  cracks. 

The  GSCM  accounts  for  the  interactions  between  cracks  in  an  average  sense  under  the  assumption  that 
each  crack  is  surrounded  by  a  region  of  intact  material.  Therefore,  in  physical  situations  when  cracks  are  close 
together,  the  strong  interactions  between  them  cannot  be  accounted  for  within  the  framework  of  the  GSCM. 
For  example,  Kachanov  [11,12]  shows  that  in  some  pre-arranged  crack  configurations,  the  crack  density  could 
become  increasingly  large,  while  having  very  little  impact  on  the  effective  moduli.  For  another  configuration 
with  a  very  low  crack  density,  the  effective  modulus  can  approach  zero.  Kachanov  [11,12]  points  out  that  the 
limitations  of  the  GSCM  stem  from  the  fact  that  this  method  does  not  permit  overlap  of  the  inclusions  con¬ 
taining  cracks.  This  does  not  preclude  modeling  of  randomly  located  cracks,  so  long  as  the  inclusions  do  not 
overlap,  yet  the  maximum  crack  density  is  limited  to  r]  —  0.3.  In  this  statement,  Kachanov  refers  to  Aboudi 
and  Benveniste  [1],  where  the  authors  use  a  circular  inclusion,  which  geometrically  limits  the  crack  density 


Fig.  1.  GSCM  model  from  Aboudi  and  Benveniste  [1]. 


Fig.  2.  Model  of  the  crack  problem  from  Santare  et  al.  [19]. 
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to  rj  =  l/n  as  discussed  above.  At  this  density,  the  crack  tips  touch  the  inclusion-matrix  interface.  In  Huang 
et  al.  [18]  and  Santare  et  al.  [19],  the  authors  use  an  elliptical  rather  than  a  circular  inclusion.  With  the  elliptical 
inclusion  this  geometric  limitation  is  eliminated  and  any  rj  value  can  be  used  in  the  calculations.  However  the 
assumption  of  no  overlapping  between  the  inclusions  is  still  present,  so  the  model  cannot  account  for  cracks 
which  nearly  touch.  Because  of  this,  there  is  still  a  practical  limitation  on  the  value  for  r\. 

In  cases  where  there  is  strong  interaction  among  multiple  cracks,  the  GSCM  is  not  useful  and  one  would 
need  other  techniques  to  calculate  the  effective  moduli.  Such  techniques  have  been  developed.  For  example, 
Greengard  and  Helsing  [20]  employ  the  fast  multipole  method  (FMM)  to  evaluate  the  stress  field  of  the  elastic 
media  with  different  inclusions;  Helsing  [21]  applies  FMM  in  calculating  the  interactions  among  large  numbers 
of  cracks;  Wang  and  Chau  [22]  use  a  boundary  element  approach  to  study  the  effective  stress  intensity  factor 
and  the  interaction  between  cracks  and  holes;  Pan  [23]  derives  the  single  domain  BEM  formulation  to  study 
the  effect  of  different  cracks  in  different  size  of  domains;  Wang  and  Feng  [24]  discuss  the  interaction  between 
multiple  rows  of  periodical  cracks;  Dong  and  Lee  [25]  apply  the  boundary  integral  equation  method  to  eval¬ 
uate  the  interaction  between  the  doubly  periodic  array  of  cracks  and  their  effect  on  the  mechanical  properties. 
These  methods  may  be  used  to  account  for  strong  crack  interactions  resulting  from  high  crack  densities  and 
presumably,  from  randomly  located  cracks.  But  the  example  solutions  in  the  referenced  literature  still  do  not 
account  for  randomly  located  cracks.  In  the  current  paper,  we  show  the  results  for  crack  density  up  to  r]  =  0.4, 
which  as  described  in  Kachanov  [11]  “can  be  considered  as  quite  high”  for  a  2-D  problem,  but  seems  a  rea¬ 
sonable  upper  limit  in  light  of  the  considerations  discussed  above. 

In  most  of  the  above  studies,  the  authors  consider  either  randomly  oriented  cracks  or  aligned  cracks.  For 
randomly  oriented  cracks,  the  effective  medium  will  be  isotropic,  since  there  is  no  preferred  crack  orientation. 
In  the  case  of  aligned  cracks,  the  effective  moduli  will  clearly  be  anisotropic,  since  the  material  behaves  differ¬ 
ently  when  the  cracks  are  open  than  when  they  are  not.  In  real  applications,  however,  the  orientation  distri¬ 
bution  of  the  cracks  will  be  neither  random  nor  aligned.  For  example,  the  cracks  shown  in  Fig.  3  have  an 
average  orientation  parallel  to  the  x-axis  and  are  distributed  within  ±0O  relative  to  the  x-axis.  In  this  case, 
crack  opening  will  influence  the  moduli  in  all  directions,  but  not  equally.  Santare  et  al.  [19]  extended  Aboudi 
and  Benveniste’s  [1]  GSCM  to  study  the  2-D  anisotropic  case  where  the  cracks  have  a  prescribed  orientation 
distribution  </>($)  (Fig.  3).  They  used  the  solution  for  a  crack,  inside  an  elliptical  inclusion  surrounded  by  an 
anisotropic  effective  medium  (Fig.  2).  Feltman  and  Santare  [26]  extended  the  Santare  et  al.  [19]  work  to  study 
the  problem  of  arbitrarily  oriented  cracks  in  a  material  that  is  originally  anisotropic  before  damage  occurs. 

In  recent  years,  the  research  in  this  field  has  continued.  Zheng  and  Du  [27]  proposed  the  effective  self-con¬ 
sistent  method  (ESCM)  and  its  simplified  explicit  version,  the  interaction  direct  derivative  method  (IDDM). 
These  methods  are  based  on  an  effective  stress/strain  field  as  in  the  MTM  mentioned  earlier.  Feng  et  al.  [28] 
proposed  a  quasi-micromechanical  method  (QMM)  to  calculate  the  overall  constitutive  relations  for  brittle 
materials  with  interacting  and  growing  microcracks.  They  combined  phenomenological  observations  with 
micromechanical  models  to  realize  some  of  the  benefits  from  both  approaches.  They  built  the  model  from 
a  micromechanical  analysis  while  using  the  orientation  domain  of  microcrack  growth  (DMG)  from  continuum 
damage  mechanics  together  with  the  crack  density  parameter,  to  characterize  the  microcrack  damage.  In 
addition  to  these  quasi-analytical  studies,  several  authors,  (see  for  e.g.,  [11,29-32]),  have  provided  numerical 
solutions  to  determine  the  effective  moduli  in  2-D  cracked  media. 

It  is  clear  that  a  microcracked  medium  will  behave  differently  in  tension  and  compression,  hence  the 
effective  anisotropic  moduli  of  such  a  medium  should  reflect  this  load-induced  anisotropy.  Horii  and 
Nemat-Nasser  [3]  discuss  the  response  of  a  medium  with  randomly  distributed  cracks  under  overall 


Fig.  3.  Cracks  distributed  within  ±0o.  0  is  the  angle  between  the  crack  and  x-axis  (label  the  x-axis  in  Fig.  2). 
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compressive  forces  using  the  SCM.  In  this  solution,  some  of  the  pre-existing  cracks  close,  some  open,  which 
lead  to  the  anisotropic  response.  They  show  the  effective  moduli  under  several  different  loading  conditions. 
However,  the  results  are  limited  to  randomly  distributed  cracks  and  the  limitations  of  the  SCM,  such  as  under 
estimating  the  moduli  and  predicting  zero  moduli,  are  still  present  in  the  results. 

The  problem  of  crack  face  overlap  has  also  been  recently  addressed  by  Wang  and  Sun  [33]  in  their  work 
using  boundary  element  methods.  In  the  present  work,  we  extend  the  GSCM  developed  in  Santare  et  al. 
[19]  and  use  computational  finite  element  methods  to  calculate  the  anisotropic  effective  moduli  for  2-D 
cracked  media  in  plane  stress,  by  considering  the  different  behaviors  of  cracks  in  tension  and  compression. 


2.  Generalized  self-consistent  method  for  anisotropic  effective  media 


To  begin  the  discussion,  consider  the  application  of  the  GSCM  to  the  determination  of  the  effective  moduli 
in  a  microcrack  damaged  solid.  For  a  solid  body  containing  discontinuities  or  voids,  the  total  strain  will  be  the 
sum  of  the  deformation  of  the  material  itself  plus  the  changes  in  the  void  volume.  For  a  specific,  external  load, 
the  average  deformation  of  the  material  is  related  to  the  average  stress  and  therefore  can  easily  be  found 
through  the  undamaged  material’s  constitutive  relations.  To  determine  the  change  in  void  volume,  one  needs 
to  know  the  initial  shape  and  distribution  of  the  voids.  In  the  current  model,  the  voids  are  assumed  to  be 
homogeneously  distributed  microcracks  and  the  strain  associated  with  them  is  therefore  related  to  the  crack 
opening  displacements.  Assuming  homogeneous  loading  and  elastic  material  response,  the  strain  relationship 
described  above  can  be  written  in  terms  of  strain  energies  ([1];  see  also  Walsh  [34]  who  derives  a  similar  expres¬ 
sion  for  isotropic  solids  using  the  reciprocal  theorem), 


2^/44  =  2^4/4  +  2V 


(2) 


In  this  expression,  S*Jki  is  the  effective  compliance  of  the  damaged  material,  and  Sijk£  is  the  compliance  of  the 
undamaged  material,  <7?.  is  the  homogeneous,  applied  Cauchy  stress  field,  and 

(3) 

is  the  traction  that  would  be  present  along  the  crack  faces  if  the  material  were  not  damaged.  The  term  in 
brackets  [ut]  is  the  crack  opening  displacement  so  that  the  integral  in  the  expression,  taken  over  the  crack  sur¬ 
face  Ck  of  each  of  the  M  cracks  in  the  volume  V,  gives  the  total  work  associated  with  crack  opening.  This 
integral  is  evaluated  to  take  into  account  the  crack  distribution,  by  considering  all  the  existing  microcrack  ori¬ 
entations,  appropriately  weighted.  To  determine  “TV”  distinct  effective  moduli,  “TV”  different  homogeneous 
stress  fields  c?.  are  applied  to  Eq.  (2),  and  the  different  forms  of  the  equation  are  solved  simultaneously  to 
determine  the  “TV”  effective  elastic  constants  (see  for  example  [19]). 


3.  Crack  face  contact 


The  crack-opening  energy  integrals  used  in  the  previous  studies  are  based  on  linear  elastic  fracture  mechan¬ 
ics  (LEFM),  which  does  not  explicitly  account  for  crack  face  contact.  Therefore,  the  strain  energy  associated 
with  compressive  loading  is  the  same  as  that  associated  with  the  tensile  loading  due  to  the  fact  that  in  LEFM 
the  crack  faces  overlap  in  compression.  In  reality,  the  crack  faces  will  open  or  close  depending  on  the  applied 
loads,  resulting  in  very  different  strain  energies  in  tension  and  compression.  Consider  a  two-dimensional  plane 
stress  orthotropic  elastic  material  with  cracks,  under  bi-axial  loads  P  and  P*  (Fig.  4a)  and  a  small  load  incre¬ 
ment  A q,  which  is  not  necessarily  a  bi-axial  load. 

We  will  discuss  the  following  three  different  cases  as  shown  in  Fig.  4b.  The  o\  and  o2  in  Fig.  4b  represent  the 
principle  stresses  in  the  2-D  plane. 

Case  1 :  P>  0,  P*  >  0  and  \P\  >  \Aq\,  \P*\  >  \Aq\ 

In  this  case,  the  pre-existing  bi-axial  tensile  loads  P  and  P*  are  much  larger  than  the  load  increment  A q. 
Because  of  this,  all  cracks  are  open  and  there  is  no  crack  face  contact.  The  results  in  Santare  et  al.  [19]  are 
valid  in  this  case. 
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Fig.  4.  (a)  A  cracked  medium  under  bi-axial  load,  (b)  Three  different  load  cases,  (c)  Three  cases  with  pure  shear  load  increment. 


Case  2:  P  <  0 ,P*  <  0  and  |P|  »  \Aq\,  |P*|  \Aq\ 

In  this  case,  the  pre-existing  bi-axial  loads  P  and  P*  are  compressive  and  much  larger  than  the  load  increment 
Aq.  Therefore,  all  cracks  remain  closed  in  mode  I  (normal).  However  the  cracks  may  still  open  in  mode  II  (shear). 

Case  3:  P  =  P*  =  0 

The  medium  is  initially  stress  free  before  the  load  increment  is  applied.  In  this  case,  the  cracks  may  either 
open  or  close  in  both  mode  I  and  II  depending  on  the  direction  of  the  load  increments  Aq. 

Fig.  4c  shows  schematically  the  cracks’  response  to  these  three  cases  under  a  pure  shear  load  increment  as 
an  example  of  the  Aq.  In  Case  1,  all  the  cracks  open;  in  Case  2,  all  the  cracks  close  in  mode  I;  and  in  Case  3, 
some  cracks  open  and  some  close. 

For  Cases  1  and  2,  four  independent  material  constants  are  required  to  fully  characterize  S*Jki  in  the  plane 
(plane  stress);  thus  four  independent  load  cases  are  needed  as  shown  schematically  in  Fig.  5  (for  Case  1)  and 
Fig.  6  (for  Case  2).  Taking  each  load  case  separately  and  assuming  that  the  undamaged  material  is  isotropic, 
the  following  equations  result  [19]: 


1  1  \  C 

——pi  A  —  —pi A  +  -  /  [ui]Pli%1  d Cfr, 

2 Eg1  2 EVl  2  {-(Jc) 

1  1  1  M  r 


(4) 

(5) 
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|  P 


Fig.  5.  Four  load  cases  for  medium  in  tension. 


|p  |p 


Fig.  6.  Four  load  cases  for  medium  in  compression. 
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(6) 

(7) 


In  these  equations,  p i,  p2 ,  q  and  h  are  the  homogeneous  load  increments  as  shown  in  Figs.  5  and  6.  E*2,  /t* 

and  JC  are  the  effective  moduli  of  the  damaged  material.  E\  and  E\  are  the  elastic  moduli  in  the  principal  direc¬ 
tions,  ii  ,  is  the  shear  modulus,  and  JC  is  the  effective  2-D  plane  stress  bulk  modulus.  These  moduli  are  related 
to  Poisson’s  ratio  through  the  equation  [19], 

K*  =  E^2  (K) 

E\  +  E*2(\  —  2v\)  ’  V  ; 

where  v\  is  one  of  the  effective  Poisson’s  ratios  of  the  damaged  material.  E,  n  and  K  are  the  undamaged  mate¬ 
rial  moduli  where  K  is  analogous  to  K*  for  the  case  of  an  isotropic  material. 

In  Case  3,  the  cracked  medium  responds  differently  under  tensile  and  compressive  loads,  hence  there  will  be 
different  Young’s  moduli,  and  bulk  moduli  under  tension  and  compression  (Fig.  7).  We  assume  the  crack  dis¬ 
tribution  is  symmetric  with  respect  to  the  x-axis.  Therefore,  when  we  shear  the  medium  relative  to  this  axis, 
half  of  the  cracks  will  open  and  half  will  close  in  mode  I  (as  depicted  in  Fig.  4c),  so  there  is  only  one  effective 
shear  modulus  /i*  in  Case  3.  Seven  independent  material  constants  are  thus  required  to  describe  this  2-D  ortho¬ 
tropic  plane  stress  model 
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In  these  equations,  p  1T,  p iC,  Pit,  Pic,  01,  hT,  hc  represent  the  homogeneous  load  increments  shown  in  Fig.  7. 
E\t  and  E*2T  are  the  effective  Young’s  moduli  under  tensile  load,  and  E\c ,  E2C  are  the  effective  Young’s  moduli 
under  compressive  load.  K*c  are  the  effective  2-D  plane  stress  bulk  moduli  under  tensile  and  compressive 
loads.  By  generalizing  the  results  of  Taya  and  Mura  [35]  and  Feltman  and  Santare  [26],  they  can  be  written  as 


K*  = 


ElTE2T 


K*  = 


^1C^2C 


■  2  if  : 
au12TJ 


c  E\c  +  E*2C(1-2v 


12C  ) 


(16) 

(17) 


where  V12T  represents  the  effective  Poisson’s  ratios  of  the  damaged  material  under  tensile  loads,  and  V*12C  is  for 
compressive  loads. 


4.  Finite  element  model 

The  results  in  Santare  et  al.  [19]  are  valid  only  for  Case  1,  as  discussed  in  the  previous  section.  In  that  case, 
the  crack  surfaces  will  never  overlap  due  to  the  applied  bi-axial  tensile  loads.  In  order  to  prevent  the  unphys¬ 
ical  crack  face  overlap  under  compression  loads  in  LEFM,  and  discuss  other  cases  mentioned  in  previous  sec¬ 
tion,  we  have  developed  a  finite  element  model  to  calculate  the  energy  associated  with  the  crack  opening 
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displacements.  Nonlinear  contact  elements  are  used  along  the  crack  surfaces,  so  that  overlap  is  eliminated,  and 
crack  face  friction  can  be  added  easily. 

The  finite  element  model  was  designed  to  mimic  as  closely  as  possible  the  analytical  boundary  value  prob¬ 
lem  used  in  the  GSCM.  This  model  is  shown  schematically  in  Fig.  8.  The  use  of  finite  elements  introduces  sev¬ 
eral  new  issues.  For  example,  the  analytical  solution  domain  used  in  GSCM  is  an  infinite  medium.  To  simulate 
this  using  an  FE  model  we  need  a  large  enough  domain  so  that  the  embedded  crack  does  not  interact  with  the 
external  boundary.  However,  the  larger  the  model,  the  more  calculation  time  needed  to  converge  to  a  set  of 
effective  moduli.  A  second  issue  is  the  interface  between  the  ellipse  and  the  surrounding  medium.  In  the  FE 
model,  the  mesh  needs  to  be  sufficiently  refined  to  simulate  the  smooth  elliptical  interface  in  the  analytical 
solution.  A  comprehensive  meshing  study  was  undertaken  and  the  numerical  results  were  compared  to  known 
analytical  solutions  [36]  to  resolve  these  domain  size  and  finite  element  mesh  refinement  issues. 

The  resulting  FE  model  is  constructed  with  2-D  plane  stress,  isoparametric  elements.  As  in  the  analytical 
solution,  the  material  inside  the  ellipse  is  isotropic  and  homogeneous  and  the  material  outside  is  a  homogeneous 
ortho  tropic  effective  material  with  a  known  principal  axis  orientation.  Note  that  in  Fig.  8,  the  crack  shown  is 
horizontal.  In  order  to  calculate  the  crack  opening  displacement  for  cracks  at  other  angles,  we  simply  modify 
the  externally  applied  loads  as  calculated  from  Mohr’s  circle.  The  external  boundary  length  of  the  FE  model  is 
L  =  500  and  the  focal  length  of  the  ellipse  equals  the  half  crack  length  c  =  30.  In  order  to  model  a  specific  crack 
density  r\,  the  crack  length  is  kept  constant,  and  the  size  of  the  ellipse  is  varied.  By  using  Eq.  (1),  together  with 
the  geometrical  relations  for  an  ellipse,  A  =  nab  and  c2  =  a2  —  b 2,  we  can  determine  a  unique  a  and  b ,  the  semi¬ 
major  and  semi-minor  axes  of  the  ellipse  for  a  given  rj.  For  the  results  that  follow,  the  FE  model  has  31,882 
elements  in  the  entire  2-D  plane  stress  model,  with  4054  elements  inside  the  ellipse  when  rj  =  0.4. 

The  effective  moduli  for  the  cracked  media  are  obtained  using  the  FE  solutions  for  the  crack  energies  in  the 
method  described  in  Section  2  and  3.  We  use  ABAQUS  6.4.4  [37]  as  the  FE  solver  and  the  program  is  written 
in  the  ABAQUS  script  language  Python.  For  the  initially  stress-free  case  (Case  3),  we  wrote  a  user  subroutine 
to  define  a  new  2-D  material  with  seven  material  constants.  On  a  PC  with  a  Pentium  4  2.0  GHz  cpu  and  1  GB 
RAM,  it  took  from  30  min  to  3  h  for  a  single  iteration  with  a  specific  r\  and  60,  depending  on  the  number  of 
unknown  material  properties  (4  for  Cases  1  and  2  and  7  for  Case  3),  whether  the  friction  existed  and  whether 
we  used  the  user  defined  material.  There  are  about  5-15  iterations  in  a  “typical”  solution,  when  we  use  the 
virgin  material  properties  as  the  initial  guess. 

To  describe  the  algorithm  of  the  FE  script,  consider  the  medium  in  tension  (Case  1)  as  an  example.  Let  us 
assume  that  the  crack  distribution  function  is  </>($)  and  the  individual  crack  angles  are  distributed  in  the  inter¬ 
val  [— 0o,6o]  relative  to  the  x-axis.  With  this,  we  can  replace  the  summations  in  Eqs.  (4)— (7)  with  the  following 
integrals  (adapted  from  [19]): 
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Fig.  8.  FEM  model  for  single  crack  inside  an  elliptical  inclusion. 


D.  Su  et  al.  /  Engineering  Fracture  Mechanics  74  (2007)  1436-1455 


1445 


Fig.  9.  Flow  chart  for  calculating  effective  moduli  of  the  medium  in  tension. 
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Fig.  9  is  a  flow  chart  of  the  algorithm  used  for  calculating  the  effective  moduli  of  the  medium  in  tension.  The 
quantities  assumed  known,  are  the  crack  density,  r\  (hence  the  elliptical  geometry  and  crack  length)  the  undam¬ 
aged  material  moduli  E,  [i  and  K  and  the  crack  orientation  distribution,  </>($),  and  range  [—do,  doi-  An  initial 
guess  is  made  for  the  four  effective  moduli  in  the  surrounding  matrix.  By  applying  these  and  the  traction 
boundary  conditions  to  the  model  as  shown  in  Fig.  5,  we  calculate  the  crack  opening  displacement  [uf]  in  each 
load  case  for  the  full  range  of  crack  orientations  considered.  These  displacements  allow  us  to  approximate  the 
double  integrals  in  Eqs.  (18)— (21).  We  can  then  use  these  four  equations  to  calculate  a  new  set  of  the  four  effec¬ 
tive  moduli  E\,  E*2,  lE  and  K*.  These  new  moduli  are  then  compared  to  the  initial  estimates  to  determine  if 
another  iteration  is  necessary  to  converge  within  a  user-defined  level  of  accuracy. 


5.  Computational  results 

As  discussed  before,  different  crack  distribution  functions  </>($)  will  result  in  different  orthotropic  effective 
moduli.  In  the  following  solutions,  we  choose  </>($)  to  be  a  constant  distribution  between  [—60,60]  as  shown  in 
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Fig.  10.  Constant  distribution  of  cracks. 

Fig.  10.  This  is  chosen  so  that  when  0O  =  90,  the  cracks  are  randomly  distributed  as  is  the  case  in  most  of  the 
literature.  Other  distributions  could  be  used  simply  by  inputting  the  proper  function  into  Eqs.  (18)— (21).  We 
assume  that  the  intact  Poisson’s  ratio  is  1/3  and  that  the  crack  tips  are  placed  at  the  foci  of  the  ellipse  in  Fig.  8. 

In  the  tables,  we  list  the  results  for  the  effective  moduli  for  all  three  load  cases  described  in  Section  3.  The 
effective  moduli  for  the  medium  in  tension  (Case  1)  are  all  within  4%  relative  error  when  compared  to  the 
results  in  Santare  et  al.  [19].  This  verifies  that  the  FE  model  developed  approximates  the  analytical  solution 
to  within  an  acceptable  level  of  accuracy  when  all  the  cracks  are  open.  In  Figs.  11  and  12,  we  compare  our 
current  results  with  several  classic  results  for  the  effective  moduli  of  cracked  media  in  2-D  plane  stress  models. 
Fig.  11  is  the  effective  Young’s  modulus  for  the  medium  with  randomly  distributed  cracks  (0O  =  9 0)  and 
Fig.  12  shows  the  effective  Young’s  modulus  in  the  direction  perpendicular  to  an  aligned  array  of  cracks 
(0o  =  O). 

As  discussed  above,  in  a  2-D  orthotropic  plane  stress  model,  four  material  moduli  (Cases  1  and  2)  or  seven 
material  moduli  (Case  3)  are  required  to  fully  define  the  material  compliance.  In  plane  stress,  the  mechanical 
properties  in  the  out-of-plane  direction  are  uncoupled  from  those  in  the  in-plane  directions  and  are  therefore 
completely  arbitrary. 

However,  in  the  literature,  many  authors  predict  the  effective  moduli  for  the  cracked  medium  using  a  2-D 
plane  strain  model.  In  plane  strain,  the  out-of-plane  compliances  affect  the  in-plane  moduli.  Yet,  many 
authors  do  not  report  the  mechanical  properties  in  the  out-of-plane  direction,  and  this  makes  it  difficult  to 
compare  our  result  to  theirs.  Alternatively,  it  is  possible  to  use  the  current  procedure,  along  with  an  assumed 
set  of  out-of-plane  properties,  to  generate  an  analogous  set  of  plane  strain  results. 


- Non-interacting  Kachanov  [11] 


Fig.  11.  Normalized  effective  Young’s  modulus  for  randomly  distributed  cracks  vs.  crack  density. 
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Fig.  12.  Normalized  effective  Young’s  modulus,  in  the  direction  perpendicular  to  the  cracks,  for  aligned  cracks  vs.  crack  density. 


The  data  in  Tables  are  plotted  in  a  series  of  contour  plots  (Figs.  13-15)  to  facilitate  visualization  and 
comparison.  In  each  contour  plot,  the  horizontal  axes  are  the  limits  of  the  orientation  distribution,  60  from  0° 
to  90°  and  the  crack  density,  r\  from  0  to  0.4.  The  vertical  axis  is  the  effective  modulus  normalized  by  the  cor¬ 
responding  undamaged  modulus.  Figs.  13  and  14  show  the  moduli  for  the  medium  in  tension  (Case  1)  and 
compression  (Case  2).  Comparing  Fig.  13  with  Fig.  14,  we  see  that  the  medium  is  generally  stiffer  under  com¬ 
pressive  loads  than  under  tensile  loads.  This  is  because  the  cracks  can  only  experience  mode  II  (shear)  opening 
under  compressive  loads  and  therefore  the  energy  dissipated  through  crack  opening  is  smaller  than  that  under 
tensile  loads  where  modes  I  and  II  are  both  possible.  As  shown  in  Fig.  14,  the  2-D  effective  bulk  modulus  is 
unchanged  under  compressive  loads.  This  is  due  to  the  fact  that  when  the  medium  is  loaded  in  bi-axial  com¬ 
pression,  all  the  cracks  will  be  closed  (in  both  normal  and  shear  modes). 

Fig.  15  shows  how  the  effective  moduli  of  the  initially  stress-free  medium  (Case  3)  vary  with  0O  and  r\.  In  this 
figure,  when  60  =  90°,  the  crack  orientations  are  random  and  therefore  the  medium  is  isotropic.  The  classical 
relationship  among  isotropic  material  moduli  is  shown  as 


E 

^  2(l+i?)’ 


(22) 


However,  the  effective  properties  in  Fig.  15  do  not  satisfy  Eq.  (22).  In  an  isotropic  medium,  where  the  tensile 
properties  are  different  from  the  compressive  properties,  the  following  relationship  can  be  easily  derived: 


EqEj 

Eq(1  +  i?t)  +Ft(1  +  i?c) 5 


(23) 


where  Ec ,  ET  are  Young’s  moduli  under  tensile  and  compressive  loads  and  vT,  vc  are  Poisson’s  ratios  in  these 
two  load  cases.  In  this  case,  the  Poisson’s  ratios  in  Eq.  (23)  can  be  calculated  from  the  2-D  bulk  moduli  and 
Young’s  moduli  through  the  following  relations: 


.  _  Ej 

T  =  2(1  —  vT)  ’ 
-  =  Ec 

C  2(1  -  vc)  ’ 


(24) 


Eqs.  (24)  and  (25)  are  the  isotropic  versions  of  Eqs.  (16)  and  (17). 


(25) 
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Table  1 

Effective  moduli  of  a  medium  in  tension 


90 

E\/E 

E\  /E 

nh/n 

k*  Ik 

do 

E\/E 

E\/E 

n\2/n 

k*  Ik 

rj  =  0.1 

rj  =  0.2 

0 

1.000 

0.553 

0.796 

0.623 

0 

1.000 

0.316 

0.633 

0.381 

10 

0.993 

0.556 

0.795 

0.623 

10 

0.985 

0.319 

0.632 

0.382 

20 

0.974 

0.566 

0.795 

0.625 

20 

0.944 

0.328 

0.630 

0.384 

30 

0.944 

0.581 

0.794 

0.627 

30 

0.883 

0.344 

0.628 

0.388 

40 

0.907 

0.601 

0.794 

0.629 

40 

0.812 

0.366 

0.626 

0.391 

50 

0.866 

0.626 

0.793 

0.631 

50 

0.739 

0.393 

0.624 

0.395 

60 

0.826 

0.653 

0.793 

0.633 

60 

0.670 

0.424 

0.624 

0.397 

70 

0.789 

0.681 

0.793 

0.634 

70 

0.611 

0.459 

0.624 

0.399 

80 

0.758 

0.708 

0.794 

0.634 

80 

0.562 

0.494 

0.625 

0.400 

90 

0.732 

0.732 

0.794 

0.634 

90 

0.526 

0.526 

0.626 

0.400 

r\  =  0.3 

rj  =  0.4 

0 

1.000 

0.182 

0.502 

0.228 

0 

1.000 

0.106 

0.398 

0.137 

10 

0.976 

0.184 

0.500 

0.229 

10 

0.965 

0.108 

0.395 

0.138 

20 

0.910 

0.191 

0.496 

0.231 

20 

0.870 

0.113 

0.388 

0.139 

30 

0.816 

0.203 

0.491 

0.234 

30 

0.740 

0.120 

0.380 

0.140 

40 

0.711 

0.220 

0.486 

0.237 

40 

0.611 

0.131 

0.372 

0.140 

50 

0.611 

0.242 

0.482 

0.239 

50 

0.484 

0.145 

0.366 

0.139 

60 

0.523 

0.268 

0.481 

0.240 

60 

0.387 

0.163 

0.364 

0.138 

70 

0.452 

0.298 

0.481 

0.240 

70 

0.316 

0.184 

0.364 

0.137 

80 

0.399 

0.330 

0.483 

0.241 

80 

0.266 

0.208 

0.366 

0.136 

90 

0.361 

0.361 

0.485 

0.241 

90 

0.234 

0.234 

0.368 

0.136 

Table  2 

Effective  moduli  of  a  medium  in  compression 

90 

E\/E 

El/E 

n\2/n 

K*/K 

e0 

E\/E 

E\/E 

n\2/n 

IC/K 

rj  =  0.1 

rj  =  0.2 

0 

1.000 

1.000 

0.797 

1.000 

0 

1.000 

1.000 

0.639 

1.000 

10 

0.994 

0.994 

0.804 

1.000 

10 

0.987 

0.987 

0.650 

1.000 

20 

0.977 

0.977 

0.823 

1.000 

20 

0.953 

0.953 

0.680 

1.000 

30 

0.955 

0.955 

0.850 

1.000 

30 

0.909 

0.909 

0.724 

1.000 

40 

0.933 

0.933 

0.878 

1.000 

40 

0.868 

0.868 

0.772 

1.000 

50 

0.917 

0.917 

0.900 

1.000 

50 

0.838 

0.838 

0.810 

1.000 

60 

0.909 

0.909 

0.911 

1.000 

60 

0.824 

0.824 

0.830 

1.000 

70 

0.910 

0.910 

0.911 

1.000 

70 

0.825 

0.825 

0.828 

1.000 

80 

0.916 

0.916 

0.902 

1.000 

80 

0.836 

0.836 

0.813 

1.000 

90 

0.924 

0.924 

0.890 

1.000 

90 

0.851 

0.851 

0.793 

1.000 

rj  =  0.3 

rj  =  0.4 

0 

1.000 

1.000 

0.516 

1.000 

0 

1.000 

1.000 

0.418 

1.000 

10 

0.980 

0.980 

0.528 

1.000 

10 

0.973 

0.973 

0.432 

1.000 

20 

0.929 

0.929 

0.564 

1.000 

20 

0.903 

0.903 

0.470 

1.000 

30 

0.864 

0.864 

0.618 

1.000 

30 

0.819 

0.819 

0.528 

1.000 

40 

0.805 

0.805 

0.677 

1.000 

40 

0.745 

0.745 

0.595 

1.000 

50 

0.764 

0.764 

0.728 

1.000 

50 

0.693 

0.693 

0.653 

1.000 

60 

0.744 

0.744 

0.754 

1.000 

60 

0.670 

0.670 

0.684 

1.000 

70 

0.745 

0.745 

0.752 

1.000 

70 

0.671 

0.671 

0.682 

1.000 

80 

0.761 

0.761 

0.731 

1.000 

80 

0.690 

0.690 

0.657 

1.000 

90 

0.782 

0.782 

0.705 

1.000 

90 

0.715 

0.715 

0.626 

1.000 

When  friction  is  present  on  the  crack  surfaces,  the  moduli  become  path  and  load  dependent.  This  is  because 
the  friction  force  is  a  non-conservative  force,  and  therefore  the  work  of  friction  depends  on  the  sequence  of 
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Table  3 


Effective  moduli  of  initially  stress-free  media 


Oo 

E'n/E 

E\c/E 

E*2t/E 

E*2c/E 

n\2/n 

K-JK 

K-JK 

do 

E\t/E 

E\c/E 

E*2t/E 

Eh/E 

hl2/n 

k;/k 

K-JK 

n  = 

0.1 

t]  =  0.2 

0 

1.000 

1.000 

0.553 

1.000 

0.796 

0.623 

1.000 

0 

1.000 

1.000 

0.316 

1.000 

0.634 

0.381 

1.000 

10 

0.993 

0.994 

0.556 

0.994 

0.799 

0.623 

1.000 

10 

0.986 

0.987 

0.319 

0.987 

0.640 

0.382 

1.000 

20 

0.975 

0.977 

0.565 

0.977 

0.809 

0.625 

1.000 

20 

0.947 

0.953 

0.328 

0.953 

0.654 

0.385 

1.000 

30 

0.946 

0.954 

0.580 

0.954 

0.821 

0.627 

1.000 

30 

0.889 

0.909 

0.344 

0.909 

0.673 

0.389 

1.000 

40 

0.910 

0.933 

0.600 

0.933 

0.834 

0.629 

1.000 

40 

0.821 

0.867 

0.366 

0.867 

0.693 

0.394 

1.000 

50 

0.871 

0.917 

0.624 

0.917 

0.844 

0.632 

1.000 

50 

0.749 

0.837 

0.393 

0.837 

0.708 

0.398 

1.000 

60 

0.831 

0.909 

0.651 

0.909 

0.849 

0.633 

1.000 

60 

0.681 

0.823 

0.424 

0.823 

0.716 

0.401 

1.000 

70 

0.793 

0.910 

0.679 

0.910 

0.849 

0.635 

1.000 

70 

0.620 

0.823 

0.459 

0.823 

0.715 

0.403 

1.000 

80 

0.760 

0.916 

0.707 

0.916 

0.845 

0.635 

1.000 

80 

0.570 

0.835 

0.496 

0.835 

0.709 

0.404 

1.000 

90 

0.733 

0.924 

0.733 

0.924 

0.840 

0.635 

1.000 

90 

0.530 

0.850 

0.530 

0.850 

0.701 

0.404 

1.000 

n  = 

0.3 

rj  =  0.4 

0 

1.000 

1.000 

0.182 

1.000 

0.505 

0.228 

1.000 

0 

1.000 

1.000 

0.107 

1.000 

0.402 

0.137 

1.000 

10 

0.978 

0.980 

0.184 

0.980 

0.511 

0.230 

1.000 

10 

0.968 

0.973 

0.108 

0.973 

0.408 

0.138 

1.000 

20 

0.916 

0.928 

0.192 

0.928 

0.526 

0.233 

1.000 

20 

0.880 

0.902 

0.114 

0.902 

0.423 

0.141 

1.000 

30 

0.828 

0.863 

0.205 

0.863 

0.547 

0.237 

1.000 

30 

0.761 

0.816 

0.122 

0.816 

0.443 

0.144 

1.000 

40 

0.728 

0.803 

0.222 

0.803 

0.569 

0.241 

1.000 

40 

0.632 

0.739 

0.134 

0.739 

0.463 

0.146 

1.000 

50 

0.630 

0.760 

0.244 

0.760 

0.586 

0.245 

1.000 

50 

0.514 

0.686 

0.149 

0.686 

0.478 

0.147 

1.000 

60 

0.542 

0.740 

0.271 

0.740 

0.595 

0.247 

1.000 

60 

0.415 

0.661 

0.168 

0.661 

0.486 

0.146 

1.000 

70 

0.469 

0.741 

0.302 

0.741 

0.594 

0.248 

1.000 

70 

0.339 

0.662 

0.191 

0.662 

0.485 

0.145 

1.000 

80 

0.412 

0.756 

0.336 

0.756 

0.587 

0.248 

1.000 

80 

0.283 

0.680 

0.217 

0.680 

0.479 

0.144 

1.000 

90 

0.370 

0.777 

0.370 

0.777 

0.578 

0.248 

1.000 

90 

0.246 

0.706 

0.245 

0.706 

0.470 

0.143 

1.000 

Table  4 

Effective  moduli  of  initially  stress-free  media  with  friction  coefficient  equal  to  0.2 

0o 

E\j/E 

Eh/E 

E*2t/E 

Eh/E 

hu  //* 

K-JK 

K-JK 

e_o_ 

E\t/E 

E\c/E 

E*2t/E 

Eh/E 

hl2/n 

K-JK 

K/K 

i]  = 

0.1 

1  = 

0.2 

0 

1.000 

1.000 

0.553 

1.000 

0.796 

0.623 

1.000 

0 

1.000 

1.000 

0.316 

1.000 

0.634 

0.381 

1.000 

10 

0.993 

0.994 

0.556 

1.000 

0.802 

0.623 

1.000 

10 

0.986 

0.987 

0.319 

1.000 

0.644 

0.382 

1.000 

20 

0.975 

0.978 

0.565 

0.994 

0.814 

0.625 

1.000 

20 

0.947 

0.955 

0.329 

0.988 

0.661 

0.385 

1.000 

30 

0.946 

0.958 

0.580 

0.978 

0.828 

0.627 

1.000 

30 

0.890 

0.915 

0.344 

0.955 

0.683 

0.390 

1.000 

40 

0.910 

0.940 

0.600 

0.959 

0.841 

0.629 

1.000 

40 

0.821 

0.880 

0.366 

0.917 

0.703 

0.394 

1.000 

50 

0.871 

0.928 

0.624 

0.943 

0.850 

0.632 

1.000 

50 

0.750 

0.857 

0.393 

0.886 

0.717 

0.399 

1.000 

60 

0.831 

0.924 

0.651 

0.934 

0.855 

0.634 

1.000 

60 

0.681 

0.850 

0.425 

0.868 

0.725 

0.402 

1.000 

70 

0.793 

0.928 

0.680 

0.932 

0.855 

0.635 

1.000 

70 

0.621 

0.856 

0.460 

0.864 

0.725 

0.404 

1.000 

80 

0.760 

0.935 

0.708 

0.936 

0.851 

0.635 

1.000 

80 

0.570 

0.869 

0.496 

0.871 

0.720 

0.405 

1.000 

90 

0.733 

0.942 

0.733 

0.942 

0.846 

0.635 

1.000 

90 

0.531 

0.883 

0.531 

0.883 

0.711 

0.405 

1.000 

n  = 

0.3 

rj  = 

0.4 

0 

1.000 

1.000 

0.182 

1.000 

0.505 

0.229 

1.000 

0 

1.000 

1.000 

0.107 

1.000 

0.413 

0.138 

1.000 

10 

0.978 

0.980 

0.184 

1.000 

0.515 

0.230 

1.000 

10 

0.968 

0.972 

0.109 

1.000 

0.418 

0.139 

1.000 

20 

0.916 

0.930 

0.192 

0.981 

0.535 

0.233 

1.000 

20 

0.881 

0.904 

0.114 

0.974 

0.436 

0.141 

1.000 

30 

0.828 

0.871 

0.205 

0.931 

0.559 

0.238 

1.000 

30 

0.762 

0.824 

0.123 

0.905 

0.458 

0.144 

1.000 

40 

0.729 

0.819 

0.222 

0.874 

0.581 

0.242 

1.000 

40 

0.634 

0.758 

0.135 

0.828 

0.479 

0.147 

1.000 

50 

0.631 

0.786 

0.245 

0.827 

0.597 

0.245 

1.000 

50 

0.515 

0.716 

0.150 

0.768 

0.493 

0.147 

1.000 

60 

0.543 

0.775 

0.272 

0.801 

0.606 

0.248 

1.000 

60 

0.417 

0.703 

0.169 

0.735 

0.501 

0.147 

1.000 

70 

0.470 

0.784 

0.303 

0.796 

0.606 

0.249 

1.000 

70 

0.340 

0.713 

0.192 

0.728 

0.501 

0.146 

1.000 

80 

0.413 

0.804 

0.337 

0.806 

0.600 

0.249 

1.000 

80 

0.285 

0.738 

0.219 

0.742 

0.495 

0.145 

1.000 

90 

0.371 

0.824 

0.371 

0.824 

0.589 

0.248 

1.000 

90 

0.247 

0.764 

0.247 

0.764 

0.486 

0.144 

1.000 

loading.  To  simplify  the  current  discussion,  here  we  will  only  consider  the  initially  stress-free  medium  with  the 
external  loads  applied  proportionally. 
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Fig.  13.  Effective  moduli  of  a  medium  in  tension. 


Fig.  14.  Effective  moduli  of  a  medium  in  compression. 


Fig.  16  shows  the  comparison  of  normalized  effective  Young’s  moduli  for  the  initially  stress-free  media  with 
and  without  friction  under  compressive  load.  The  crack  density  rj  =  0.4,  and  in  the  friction  case,  the  friction 
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coefficient  is  equal  to  0.2.  We  can  see  that  the  friction  affects  to  the  Young’s  moduli  under  compressive  loads. 
When  r\  =  0.4,  E\c  and  E\c  increase  about  8%  when  the  friction  exists.  Also,  we  see  that  there  is  a  clear  “flat 
zone”  (E*2C/E  =  1)  for  E\ c  when  60<  11°.  This  is  because  the  friction  prevents  the  cracks  from  opening  in  any 
mode  under  this  situation.  The  maximum  extent  of  the  “flat  zone”  is  determined  by  the  friction  coefficient. 


Schematic  cracked  Compression  Tension 

media 

Normalized  slowness  in  y  Normalized  slowness  in  y 
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Fig.  16.  Normalized  E*]ClE*2C  of  initially  stress-free  media  with  and  without  friction,  rj  =  0.4  and  the  friction  coefficient  equal  to  0.2. 


Fig.  17.  Normalized  slowness  of  cracked  media  under  tension  and  compression.  Crack  density  rj  =  0.4.  Outside  lines:  quasi- transverse 
wave;  inside  lines:  quasi-longitudinal  wave.  All  the  lines  are  normalized  by  the  slowness  of  transverse  wave  in  virgin  material. 
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6.  Wave  propagation  in  microcrack  damaged  media 


When  applying  the  predicted  effective  moduli  to  the  wave  equation,  we  can  predict  the  elastic  wave  speed 
within  the  cracked  media.  To  find  out  the  wave  speed  of  a  plane  harmonic  wave  in  the  L  direction,  we  need  to 
solve  the  following  eigenvalue  equation  [38]: 

detj^  —  Sijpv2}  =  0,  (26) 

in  which,  A  =  LTCL  is  the  acoustic  tensor,  C  is  the  elastic  tensor  and  L  is  a  3  x  2  direction  cosine  matrix. 


[L] 


l\  0 
0  l2 
h  h 


(27) 


The  slowness  is  the  inverse  of  wave  speed 

S,r’(l)  =  ^-  (28) 

Fig.  17  shows  the  normalized  slowness  (inverse  of  the  wave  speed)  in  the  media  with  cracks  in  tension  and 
compression.  The  outside  lines  depict  curves  of  maximum  slowness  (quasi-transverse  wave)  and  the  inside  lines 
depict  curves  of  minimum  slowness  (quasi-longitudinal  wave).  The  figures  on  the  first  line  are  the  media  under 
tensile  loads  (Case  1  in  Fig.  4b)  and  the  figures  on  the  second  line  are  the  media  under  compressive  loads  (Case 
2  in  Fig.  4b).  Comparing  the  figures  in  the  same  column,  we  can  see  that  for  identical  crack  distributions  and 
densities,  the  wave  speed  (slowness)  profiles  are  very  different  under  different  loads. 


7.  Conclusions 


In  this  study,  a  finite  element  model  was  used  to  solve  the  analytical  crack-inclusion  anisotropic  boundary 
value  problem  used  in  the  GSCM.  The  analytical  results  of  Santare  et  al.  [19]  were  used  to  validate  the  finite 
element  model  within  a  4%  relative  error  (Case  1),  which  confirms  the  accuracy  of  the  model.  Using  the  same 
basic  FE  model,  contact  elements  were  introduced  to  eliminate  the  crack  surface  overlap,  which  occurs  under 
compression  in  linear  elastic  fracture  mechanics.  With  this  model,  we  calculated  effective  moduli  for  a  dam¬ 
aged  medium  under  tensile,  compressive,  and  initially  stress-free  conditions.  Base  on  the  effective  moduli, 
we  also  predict  the  elastic  wave  slowness  within  the  2-D  cracked  media. 

Generally  speaking,  the  results  show  that  the  existence  of  microcracks  decreases  the  effective  stiffness  of  the 
medium  and  the  effective  moduli  further  decrease  as  the  microcrack  density  is  increased.  The  crack  orientation 
distribution  also  has  a  significant  effect  on  the  anisotropy  of  the  effective  moduli  which  can  vary  from  highly 
anisotropic  when  the  cracks  are  unidirectional  (0O  =  0°)  to  isotropic  when  they  are  randomly  distributed 
(Go  =  90°).  Furthermore,  we  see  that  the  damaged  medium  responds  quite  differently  depending  on  whether 
the  external  loads  are  tension,  compression  or  initially  zero.  And  we  see  that  the  medium  is  much  stiffer  when 
the  overall  external  loads  are  compressive. 

We  have  described  how,  in  an  initially  stress-free  medium  under  2-D  plane  stress  loading,  there  are  seven 
independent  orthotropic  effective  material  properties  for  the  microcracked  medium  (Case  3);  for  randomly  dis¬ 
tributed  microcracks,  the  number  of  independent  moduli  are  reduced  to  four,  and  satisfy  Eq.  (23).  We  have 
also  shown  how  friction  can  be  introduced  into  the  model,  and  that  in  general,  friction  makes  the  medium 
stiffer  particularly  under  compression  loads.  From  the  literature  [3]  we  know  that  the  moduli  become  path 
and  load  dependent  when  crack  face  friction  is  included  in  the  analysis,  and  to  adequately  account  for  these 
effects,  further  study  is  required. 

Finally,  we  plot  wave  slowness  profiles  for  microcrack  damaged  media  and  show  that  for  media  with  iden¬ 
tical  microcrack  densities  and  distributions  the  wave  speeds  are  very  different  under  tensile  and  compressive 
loads,  because  of  the  different  effective  moduli  for  these  two  load  cases. 
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Direct  numerical  simulations  of  waves  traveling  through  microcrack-damaged  media  are 
conducted  and  the  results  are  compared  to  effective  medium  calculations  to  determine 
the  applicability  of  the  latter  for  studying  wave  propagation.  Both  tensile  and  compressive 
waves  and  various  angular  distributions  of  randomly-located  cracks  are  considered.  The 
relationships  between  the  input  wavelength  and  the  output  wave  speed  and  output  signal 
strength  are  studied.  The  numerical  simulations  show  that  the  wave  speed  is  nearly  con¬ 
stant  when  1  /ka  >  60  for  tensile  waves  and  1  /ka  >10  for  compressive  waves,  where  k  is 
the  wave  number  and  a  is  the  average  half-crack  length.  The  direct  simulations  also  show 
that  when  the  input  wavelength  is  much  longer  than  the  crack  length,  1  / ka  >  60,  the  wave 
can  pass  through  the  damaged  medium  relatively  unattenuated.  On  the  other  hand,  when 
the  input  wavelength  is  shorter  than  a  “cut  off’  wave  length,  the  output  wave  magnitude 
decreases  linearly  with  the  input  wavelength.  The  effective  medium  wave  speed  and  mag¬ 
nitude  calculations  are  not  dependent  on  the  input  wavelength  and  therefore  the  results 
correspond  well  with  the  numerical  simulations  for  large  1  / ka .  This  suggests  a  minimum 
wavelength  for  which  the  homogenized  methods  can  be  used  for  studying  these  problems. 

©  2008  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

For  many  brittle  materials,  the  most  common  damage  mechanism  is  microcracking.  As  microcracks  develop,  they  change 
the  local  mechanical  response  to  input  loads.  One  way,  commonly  used  to  evaluate  the  quasi-static  response  of  the  damaged 
medium  is  by  directly  simulating  an  actual  array  of  cracks.  Solutions  of  this  type  include  the  fast  multipole  method  (FMM, 
[1,2]),  boundary  element  method  (BEM  [3,4])  and  many  other  numerical  solutions,  such  as  described  in  [5-9].  The  major 
advantage  of  these  methods  is  that  they  take  into  account  the  actual  crack-crack  interactions.  There  is  no  theoretical  lim¬ 
itation  to  the  crack  density  or  crack  distribution  and  these  direct  methods  can  give  highly  accurate  results  if  the  model  is 
an  accurate  representation  of  the  damaged  medium.  However,  the  detailed  calculation  of  the  interaction  between  the  cracks 
is  also  the  major  disadvantage  of  the  direct  methods.  As  the  crack  density  increases,  the  modeling  time  and  computational 
demands  can  increase  dramatically. 

To  address  this  increase  in  computational  demand,  effective  medium  methods  can  be  used.  In  these  methods,  one  models 
the  quasi-static  response  of  a  medium  with  diffuse  microcrack  damage  by  replacing  the  damaged  material  with  an  effective 
medium,  possessing  the  same  local  mechanical  properties.  Many  such  effective  medium  models  have  been  developed,  such 
as  the  self  consistent  method  (SCM)  [10-12],  Generalized  self  consistent  method  (GSCM)  [13-17],  Mori-Tanaka  method 
(MTM)  [18,19],  and  the  differential  scheme  method  (DSM)  [20].  For  a  comprehensive  review  of  the  literature  in  effective 
mechanical  properties  prediction,  readers  can  refer  to  [6,21]. 
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While  the  above-cited  literature  primarily  discusses  the  quasi-static  response  of  microcrack-damaged  media,  several 
researchers  have  addressed  elastic  wave  propagation  behavior  in  these  materials.  In  1995  Sayers  and  Kachanov  [22]  used 
the  effective  moduli  of  a  cracked  medium  to  study  microcrack-induced  elastic  wave  anisotropy.  In  that  paper,  the  authors 
used  second-rank  and  fourth-rank  crack  density  tensors  to  evaluate  the  effective  elastic  moduli.  They  calculated  the  effective 
moduli  for  the  (randomly-oriented)  isotropic  distribution  of  cracks  and  the  perfectly  aligned  cracks.  For  the  isotropic  cracked 
medium,  they  also  back-calculated  the  crack  density  from  experimental  data  of  the  ultrasonic  wave  speed.  Following  Sayers 
and  Kachanov’s  results,  Schubnel  and  Gueguen  [23,24]  used  micromechanical  and  statistical  physics  to  evaluate  the  elastic 
wave  velocities  and  permeability  of  cracked  rocks.  They  calculated  the  velocity  for  high  and  low  frequency  waves  for  both 
wet  and  dry  cracks  for  an  aligned  array  of  cracks  and  an  array  composed  of  ±60°  cracks  under  an  external  confining  pressure. 
Maurel  et  al.  [25]  studied  the  elastic  wave  propagation  through  a  random  array  of  dislocations.  Markov  [26]  used  the  Frenkel- 
Biot  theory  to  study  the  longitudinal  harmonic  wave  propagation  in  an  isotropic  porous  matrix  with  inclusions.  Levin  and 
Markov  [27]  used  an  effective  medium  approximation  (EMA),  based  on  the  self-consistent  method,  to  determine  the  effective 
elastic  moduli  and  elastic  wave  propagation  velocities  in  a  transversely  isotropic  solid  containing  aligned  inclusions. 

Yet,  with  all  this  literature,  there  are  several  questions  which  remain.  For  example,  under  what  conditions  can  we  model 
an  actual  damaged  material  with  an  equivalent  effective  medium?  How  accurate  is  the  effective  medium  representation? 
What  factors  affect  the  accuracy? 

Several  researchers  have  performed  studies  to  address  some  of  these  topics.  Zhang  and  Gross  [28]  studied  the  propagation 
of  an  elastic  wave  through  a  medium  with  randomly  distributed  cracks.  The  cracks  are  treated  as  finite  length  lines  with  dis¬ 
placement  discontinuity  and  the  crack  surfaces  are  assumed  to  be  frictionless.  They  studied  the  scattering  function  for  a  sin¬ 
gle  crack  and  expanded  this  using  a  numerical  approach  to  study  the  response  of  multiple  cracks.  They  used  both  the  theory 
of  Foldy  and  the  causal  approach  based  on  K-K  relations  to  compute  the  attenuation  coefficient  and  the  effective  wave  veloc¬ 
ity.  Littles  et  al.  [29]  performed  an  experimental  and  theoretical  investigation  to  study  how  longitudinal  waves  are  scattered 
by  a  distribution  of  cracks.  The  results  show  that  the  transmission  coefficients  are  a  function  of  the  incident  wave  number, 
the  crack  size  and  crack  density.  Kawahara  and  Yamashita  [30]  studied  the  scattering  of  P,  SV  and  SH  waves  by  a  zonal  dis¬ 
tribution  of  uniaxial  cracks  using  a  theoretical  analysis.  They  showed  that  the  attenuation  coefficient  peaks  around  ka  «  1, 
the  phase  velocity  is  almost  independent  of  k  in  the  range  of  ka  <  1  and  increases  monotonically  when  k  is  in  the  range  of 
ka  >  1.  They  also  found  that  the  effect  of  crack-face  friction  is  to  shift  the  peak  of  the  attenuation  coefficient  toward  the  lower 
wave  number  range.  Kelner  et  al.  [31]  used  the  boundary  element  method  to  simulate  P  wave  propagation  in  media  with 
uniaxial  cracks.  The  authors  conclude  that  when  the  wavelength  of  the  incident  wave  is  close  to,  or  shorter  than,  the  crack 
length,  the  scattered  waves  are  efficiently  excited.  They  also  observed  that  the  attenuation  factor  of  the  direct  P  wave  peaks 
at  around  ka  =  2,  where  k  is  the  incident  wave  number  and  a  is  the  crack  length,  and  decreases  proportionally  with  (ka)-1  in 
the  high  wave  number  range. 

Most  of  the  above  papers  focus  on  low  crack  densities  and  high  frequency  input  waves,  a  combination  which  corresponds 
to  high  attenuation.  Also,  the  crack  orientations  in  most  of  the  papers  are  limited  to  aligned  or  randomly-oriented  distribu¬ 
tions.  Additionally,  in  most  of  the  literature  cited,  the  authors  assume  the  crack  surfaces  are  stress  free,  which  means  the 
cracks  are  opened  under  all  loads.  In  our  previous  paper  [16],  we  use  the  generalized  self  consistent  method  (GSCM)  to  pre¬ 
dict  the  anisotropic  effective  moduli  of  a  medium  containing  an  arbitrarily-oriented  distribution  of  cracks.  In  that  paper,  we 
show  that  because  the  cracks  open  or  close  under  different  external  loads,  the  effective  moduli  under  different  loading  con¬ 
ditions  are  quite  different.  In  that  paper,  we  discussed  three  different  loading  conditions:  (1 )  overall  applied  tensile  loads;  (2) 
overall  compressive  loads;  (3)  initially  stress  free  media.  The  effective  moduli  under  these  three  different  loading  conditions 
are  evaluated  taking  into  account  the  effects  of  crack-face  contact.  In  that  paper  we  also  mentioned  that  tensile  and  com¬ 
pressive  wave  propagation  in  the  same  medium  behave  very  differently.  In  this  paper,  we  study  under  what  conditions 
one  may  use  the  calculated  effective  moduli  to  simulate  wave  propagation  in  the  microcrack  damaged  medium.  We  conduct 
direct  numerical  simulations  to  relate  input  wavelength  to  propagating  wave  speed  and  output  signal  strength,  and  compare 
these  results  to  those  obtained  via  the  effective  medium  models.  Once  these  comparisons  are  established,  we  can  establish 
guidelines  for  the  use  of  effective  medium  calculations  to  analyze  ultrasound  data  for  determining  the  internal  state  of 
microcrack  damage. 


2.  FEM  model 


In  the  following  numerical  experiments,  we  use  the  definition  of  crack  density  and  crack  orientation  distribution  function 
cj){6)  described  in  [15,16].  We  define  the  crack  density  as 


rj  = 


Mi i2 
A 


(1) 


where  M  is  the  number  of  cracks  per  unit  area  A,  and  a  is  the  average  half-crack  length.  We  assume  the  crack  orientation  can 
be  described  by  an  orientation  distribution  function  where  9  is  the  angle  of  the  individual  crack  with  respect  to  the 
average  crack  orientation,  assigned  to  be  the  x-axis.  In  the  following  analysis,  we  assume  that  the  cracks  are  evenly 
distributed  between  the  angles  -90  and  +60,  as  shown  in  Fig.  1.  A  more  detailed  description  of  the  crack  distribution  function 
</>($),  and  crack  density  function  rj ,  can  be  found  in  [16]. 
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crack 


Fig.  1.  Even  distribution  of  crack  orientation. 


To  simulate  uni-directional  plane  wave  propagation,  we  assume  an  infinite  half  plate,  and  apply  a  time-dependent  load 
uniformly  along  the  top  surface  (Fig.  2a).  Because  the  plate  is  infinitely  extended  in  the  x-direction,  we  can  pick  a  strip  from 
anywhere  in  the  half  plane  as  a  representative  sample  (Fig.  2b)  by  applying  symmetry  boundary  conditions  to  the  left  and 
right  edges  of  the  strip.  The  location  of  each  crack  within  the  strip  was  generated  by  a  random  number  generator  as  was  the 
crack  angle  6  chosen  to  be  in  the  interval  +60  and  -0O.  If  the  randomly  generated  location  of  a  crack  caused  it  to  overlap  with 
an  already  existing  crack,  that  crack  was  eliminated  from  the  model  and  the  process  was  continued  until  the  desired  crack 
density  was  reached.  Cracks  that  contacted  the  left  or  right  boundaries  were  kept  in  the  model,  but  because  of  the  symmetry 
conditions  imposed  on  these  boundaries,  their  lengths  could  be  quite  different  from  the  average.  Building  the  model  in  this 
manner,  the  array  is  not  completely  random,  but  represents  a  relatively  homogeneous  approximation  to  randomly  distrib¬ 
uted  cracks.  Contact  elements  are  placed  along  the  cracks  to  prevent  crack-face  interpenetration.  Fig.  2b  is  a  representation 
of  the  models  used  in  the  numerical  experiments.  A  number  of  trial  models  were  run  to  determine  an  appropriate  length  and 
width  for  the  strips  and  the  chosen  dimensions  represent  a  compromise  between  a  model  large  enough  to  eliminate  edge 
effects  and  one  small  enough  to  allow  for  efficient  computations.  All  lengths  are  normalized  with  respect  to  the  average  crack 
length  (2a)  and  the  resulting  width  of  each  model  is  4,  and  the  length  is  80.  In  all  the  results  that  follow,  the  crack  density 
parameter  is  chosen  to  be  rj  =  0.4.  (Most  of  the  previously  reviewed  literature  studies  crack  densities  in  the  range  of  0.1 -0.2 
when  defined  by  the  term  rj,  in  Eq.  (1))  Larger  densities  proved  difficult  to  model,  due  to  interference  between  individual 
cracks  and  smaller  densities  show  the  same  overall  trends,  but  to  a  lesser  extent.  Since  the  cracks  are  randomly  located 
and  oriented,  we  constructed  four  or  five  models  for  each  orientation  distribution  to  show  the  variability  introduced  by 
the  randomness. 

In  the  simulation,  we  use  ABAQUS  6.4.4  as  the  finite  element  solver.  In  Figs.  3-6,  we  show  the  finite  element  models  for 
samples  1-14.  As  shown  in  Fig.  3d,  each  crack  is  represented  by  two  crack  surfaces  with  zero  distance  between  them.  There¬ 
fore,  in  the  stress  free  state,  all  the  cracks  are  closed. 

In  order  to  study  the  relationship  between  the  input  wavelength  and  the  output  signal,  we  apply  a  single  frequency  time- 
dependent  load,  uniformly  along  the  top  of  the  strip.  Fig.  7  shows  a  sample  input  load  with  the  function 


P(t)=A 


T  -  cos  (cot) 
2 


for 


_  cot 

0<=-<l, 

2n 


(2) 


where,  P(t)  is  the  time-dependent  load,  A  is  the  maximum  load  amplitude,  co  is  the  circular  frequency  of  the  input  and  t  is  the 
time.  We  can  apply  positive  or  negative  pressure  in  order  to  generate  a  compressive  or  tensile  wave  inside  the  sample. 

In  the  numerical  experiments,  we  use  the  dimensionless  number  1  / ka  as  the  effective  wave  length  parameter,  where  a  is 
the  average  half-crack  length,  k  =  <x>/cp,  co  is  the  input  wave  frequency  and  cpis  the  plane-stress  longitudinal  wave  velocity  in 
the  un-damaged  medium  [32],  which  is  defined  as 


Fig.  2.  Model  schematics,  (a)  Time-dependent  load  is  applied  along  the  top  boundary  of  the  half  plane,  (b)  A  representative  strip  from  the  plane  is  used  in 
the  finite  element  simulation. 


Sample  04  Sample  05 


Fig.  4.  The  top  part  of  finite  element  model  of  samples  01-05  (0O  =  0°). 


CD  = 


(1  -  v2)p' 


(3) 


This  effective  wave  number  1  Ika,  therefore  represents  the  ratio  between  the  input  wave  length  in  the  undamaged  media  and 
the  average  half  crack  length. 
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Sample  06 


Sample  07 


Fig.  5.  The  top  part  of  finite  element  model  of  samples  06-09  (0O  =  30°). 


Sample  10  Sample  1 1  Sample  12 


Sample  13  Sample  14 


Fig.  6.  The  top  part  of  finite  element  model  of  samples  10-14  (0O  =  60°). 
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COt/(27l) 


Fig.  7.  Sample  input  load.  The  function  of  the  input  load  is  P  =  A 1  c°s(o}t)  for  0  ^  ^  1. 


3.  Wave  speed  calculation 

In  this  section  we  determine  the  wave  speeds  in  the  direct  numerical  simulations  and  compare  these  to  the  speeds  in  the 
effective  media  models.  The  properties  of  the  effective  media,  and  the  generalized-self-consistent  method  used  to  calculate 
them,  are  shown  in  [16].  All  the  wave  speeds  are  normalized  by  the  plane-stress  wave  speed  in  the  intact  material  cp,  shown 
in  Eq.  (3).  We  trigger  the  wave  from  the  top  of  the  models  with  time  dependent  pressure  loads  and  trace  the  stress  compo¬ 
nent  Gy.  We  chose  1%  of  the  maximum  input  as  the  threshold  value  for  the  stress  at  the  wave  front.  This  means  that  we  as¬ 
sume  that  the  wave  arrives  at  a  particular  cross  section  when  the  average  stress  oy  in  that  cross  section  reaches  1%  of  the 
maximum  input  stress. 

In  Figs.  8-10,  we  show  the  wave  speed  in  each  of  the  direct  numerical  simulation  models  as  a  function  of  the  input  wave 
length  number  1/ka,  as  well  as  the  effective  medium  predictions  and  an  analytical  solution.  The  analytical  solution  is  the 
uniaxial  longitudinal  wave  speed,  which  can  be  derived  from  [33],  for  an  anisotropic  plane  stress  medium  calculated  from 

cp2  =  aL - (4) 

y(l—  V12V21  )p 

where  E2,  vi2,  v2i  are  the  corresponding  effective  Young’s  moduli  in  y  direction,  and  the  effective  Poisson’s  ratios,  which  are 
shown  in  [16]  and  p  is  the  corresponding  mass  density.  Fig.  8  shows  results  for  the  five  different  models  with  an  aligned 
array  of  cracks  90  =  0°.  Fig.  9  contains  the  models  with  crack  orientations  evenly  distributed  between  the  values  90  =  ±30° 
and  Fig.  10,  the  models  with  angles  between  90  =  ±60°.  The  left  column  in  each  figure  shows  the  tensile  wave  speed,  while 
the  right  column  shows  the  speed  of  compressive  waves. 
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Fig.  8.  Normalized  wave  speed  vs.  1  / ka  for  90  =  0°.  Left  :  speed  of  tensile  waves.  Right  :  speed  of  compressive  wave. 
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Fig.  9.  Normalized  wave  speed  vs.  1/ka  for  90  =  30°.  Left  figure:  speed  of  tensile  waves.  Right  figure:  speed  of  compressive  wave. 


<90  =  60° 


Fig.  10.  Normalized  wave  speed  vs.  l/fca  for  0O  =  60°.  Left  :  speed  of  tensile  waves.  Right  :  speed  of  compressive  wave. 


From  the  results,  we  can  see  that  in  all  cases,  the  compressive  waves  travel  much  faster  than  the  tensile  waves.  This  is 
because  under  compression,  the  cracks  can  only  open  in  mode  II  (shear  opening)  (As  we  discuss  in  Section  2,  in  the 
stress-free  state,  the  distance  between  the  cracks  surfaces  is  zero.  Therefore,  the  crack  surfaces  cannot  move  closer  to  each 
other  under  compression.)  while  in  tension,  the  cracks  can  open  in  both  mode  I  (normal  opening)  and  mode  II.  Therefore  for  a 
given  medium,  the  material  is  much  less  stiff  when  tensile  waves  pass  through,  than  when  compressive  waves  pass  through. 
For  a  detailed  discussion  regarding  the  different  effective  material  properties  under  different  loads,  the  reader  is  referred  to 
[16].  The  figures  also  show  that  the  effective  medium  models  correspond  well  with  the  numerical  experiments  in  the  cases 
considered.  For  the  case  of  tensile  wave  speed  when  60  =  0°,  the  analytical  equations  and  the  effective  media  model  give 
higher  wave  speed  predictions  than  the  direct  numerical  simulations.  This  is  because  for  tensile  waves,  when  90  =  0°,  the 
cracks  in  the  numerical  simulations  cause  significant  energy  dispersion  as  the  wave  travels  along  the  strip.  Therefore,  for 
the  same  input  magnitude,  it  takes  longer  for  the  average  stress  across  the  strip  to  reach  the  1%  threshold  value,  resulting 
in  a  slightly  lower  calculated  speed  for  the  direct  simulations. 

Kawahara  and  Yamashita  [30]  found  analytical  solutions  for  a  similar  problem  containing  randomly  distributed,  aligned 
cracks  under  confining  pressure.  They  found  that  when  the  crack  surfaces  are  friction  free,  the  phase  velocity  is  nearly  inde¬ 
pendent  of  k  when  1/ka  >  1  and  increases  monotonically  with  k  when  1/ka  <  1.  From  the  right  column  of  Figs.  8-10  we  can 
see  that  the  wave  speed  in  our  numerical  experiments  follows  this  same  trend  and  are  completely  consistent  with  the  ana- 
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1/  ka  =  8.24  1/  ka  =  3.53  l/ka  =  0J06 


Fig.  11.  Tensile  waves  with  different  effective  wave  number  traveling  in  sample  05,  rj  =  0.4,  0o  =  0°.  Normalized  oy  is  the  average  stress  ay  in  the 
cross  session  normalized  by  the  maximum  average  stress  ay  at  distance  equal  to  zero  (the  top  boundary).  The  distance  from  the  top  is  shown  in 


M  ka-  54.9 
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Fig.  12.  Compressive  waves  with  different  effective  wave  number  traveling  in  sample  05,  r\  =  0.4, 0o  =  0°.  Normalized  oy  is  the  average  stress  ay  in  the  specific 
cross  session  normalized  by  the  maximum  average  stress  oy  at  distance  equal  to  zero  (the  top  boundary).  The  distance  from  the  top  is  shown  in  Fig.  2b. 
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lytical  results  in  Kawahara  and  Yamasita’s  paper  for  cracks  under  compression.  At  the  same  time,  from  the  left  column  in 
Figs.  8-10  we  see  that  the  tensile  wave  speed  is  nearly  constant  when  1/ka  >60  and  increases  monotonically  when  1/ 
ka  <  60.  Therefore,  we  can  see  that  the  Kawahara  and  Yamasita’s  general  conclusion  still  applies  to  tensile  wave  propagation 
in  our  current  study,  although  the  cut  off  ka  (or  1/ka)  value  is  much  different  for  tensile  waves  than  compressive  waves. 

4.  Signal  strength 

In  the  direct  numerical  simulations,  we  applied  a  uniform,  time  varying  load  along  the  top  boundary  of  the  finite  element 
models,  as  shown  in  Fig.  2b.  For  a  specific  time  and  a  specific  cross  section  in  the  model,  we  record  the  average  stress  com¬ 
ponent  Gy  through  this  cross  section  (Fig.  2b)  then  normalize  this  value  by  the  maximum  average  stress  oy  calculated  at  the 
top  boundary  of  the  model.  (Recall  that  the  crack  length  2 a  equals  1,  except  for  the  cracks  that  touch  the  boundary.)  The 
properties  of  the  intact  material  are  E//i  =  8/3,  v  =  1  / 3,  where  E  is  Young’s  modulus,  /i  is  shear  modulus  and  v  is  Poisson’s  ratio. 
The  normalized  time  is  defined  as  t  =  =  tcp/(2a). 

Fig.  1 1  shows  the  results  for  the  tensile  waves  with  different  effective  wave  numbers  traveling  in  sample  05.  We  can  see 
that  for  tensile  waves,  when  the  input  wavelength  is  much  larger  than  the  crack  length,  (1  /ka  >  10),  the  input  signal  passes 
through  the  damaged  medium  relatively  unattenuated  (the  first  row  of  the  figures).  When  the  input  wavelength  is  smaller, 
(1  /ka  =  8.24  and  3.53)  the  cracks  show  a  larger  “blocking”  effect  and  the  transmitted  wave  shows  smaller  stress  amplitude. 
When  the  input  wavelength  is  about  the  same  as  the  crack  length,  (1  /ka  =  0.706)  the  signal  is  almost  completed  blocked  by 
the  damage. 
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Fig.  13.  Output  signal  strength  vs.  input  effective  wave  number  (1  /ka)  for  the  models  with  rj  =  0.4,  0o  =  0°  and  2 a  =  1.  Left  column:  Tensile  waves  output. 
Right  column:  Compressive  waves  output.  First  row:  The  normalized  peak  averaged  stress  measured  at  distance  equal  to  10.  Second  row:  The  normalized 
peak  averaged  stress  measured  at  distance  equal  to  20. 
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Fig.  12  shows  the  compressive  wave  propagation  through  the  same  model  sample  05.  Comparing  Fig.  11  with  Fig.  12,  we 
can  see  that  the  compressive  wave  signal  does  not  decrease  appreciably,  even  when  the  input  wavelength  (1/fca)  is  relatively 
small.  This  is  because  when  tensile  waves  pass  through  the  sample,  the  cracks  are  opened  and  parts  of  the  waves  are  re¬ 
flected  by  the  open  crack  surfaces,  causing  the  output  signal  to  decrease.  Yet,  when  compressive  waves  travel  in  the  model, 
all  the  cracks  remain  closed  and  therefore  have  almost  no  effect  on  the  wave  propagation  in  this  situation. 

5.  Output  signal  strength 

5.2.  Uniform  crack  length 

In  this  section,  we  show  how  the  peak  stress  values  at  different  cross  sections  change  with  varying  input  wavelength.  In 
Figs.  13-15,  AO,  A10  and  A20  represent  the  peak  average  stress  oy  at  cross  sections  with  normalized  distances  of  0, 10  and  20, 
respectively  from  the  top  of  the  model.  Fig.  13  shows  results  for  the  models  with  an  aligned  array  of  cracks  90  =  0°  and  Figs. 
14  and  15  contain  the  models  with  cracks  oriented  between  the  values  90  =  ±  30°  and  90  =  ±  60°  respectively.  As  before,  the 
left  column  in  each  figure  shows  the  results  for  the  tensile  waves,  while  the  right  column  shows  the  results  for  the  compres¬ 
sive  waves. 

In  Figs.  13-15,  the  individual  lines  represent  the  different  model  samples.  In  all  of  these  models  rj  =  0.4  and  2a  =  1  (except 
for  the  cracks  that  touch  the  boundaries).  All  of  these  models  show  the  same  phenomenon  observed  in  the  previous  section 


0o=3O° 


Fig.  14.  Output  signal  strength  vs.  input  wave  length  number  (1  jka)  for  the  models  with  r\  =  0.4,  0o  =  30°  and  2a  =  1.  Left  column:  Tensile  waves  output. 
Right  column:  Compressive  waves  output.  First  row:  The  normalized  peak  averaged  stress  measured  at  distance  equal  to  10.  Second  row:  The  normalized 
peak  averaged  stress  measured  at  distance  equal  to  20. 
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0,  =  60° 


Fig.  15.  Output  signal  strength  vs.  input  wave  length  number  (1/fca)  for  the  models  with  rj  =  0.4,  0o  =  60°  and  2a  =  1.  Left  column:  Tensile  waves  output. 
Right  column:  Compressive  waves  output.  First  row:  The  normalized  peak  averaged  stress  measured  at  distance  equal  to  10.  Second  row:  The  normalized 
peak  averaged  stress  measured  at  distance  equal  to  20. 


when  we  discussed  the  results  of  sample  05.  What  these  figures  bring  out,  is  the  variability  from  the  different  arrays  of  cracks 
(different  experimental  samples)  as  well  as  more  clearly  delineating  the  signal  strength-wavelength  relationship. 

For  tensile  waves,  when  the  input  wave  length  is  much  larger  than  the  crack  length  (l/ka  >  20),  the  input  signal  passes 
through  the  microcrack  damaged  medium  relatively  easily.  The  output  signal  strength  remains  nearly  constant  for  all  input 
wavelengths  above  (1/ka  >  20),  and  the  output  signal  strengths  are  usually  more  than  80%  of  the  input  signal  strength.  On  the 
other  hand,  the  cracks  show  larger  “blocking”  effect  when  the  input  wave  length  decreases.  In  the  area  where  1  / ka  <15,  the 
output  signal  strength  decreases  almost  linearly  with  the  input  \\ka  decreases.  Therefore,  we  can  think  of  \\ka  =  15  as  a  “cut 
off’  wave  length  below  which  tensile  waves  will  not  travel  through  a  microcrack  damaged  medium.  When  the  input  wave¬ 
length  is  close  to  or  smaller  than  the  crack  length  (1  /ka  «  1),  almost  no  signal  can  pass  through  the  medium. 

For  compressive  waves,  we  can  see  the  similar  phenomenon  but  there  are  two  differences  that  distinguish  the  compres¬ 
sive  wave  output  from  the  tensile  wave  output.  First,  when  the  1  / ka  number  is  large,  the  output  signal  of  the  compressive 
waves  is  much  more  uniform  than  that  of  the  tensile  waves.  Secondly,  the  “cut  off’  wave  length  for  the  compressive  waves  is 
much  smaller  than  that  for  the  tensile  waves.  In  this  study  we  see  that  for  compressive  waves,  the  output  signal  decreases 
only  when  1  / ka  <  5. 

These  results  suggest  that  there  is  a  “safe  range”  for  modeling  microcrack-damaged  media  with  homogenized  effective 
medium  approximations.  When  the  input  wavelength  is  larger  than  the  “cut  off’  wavelength,  the  signal  can  travel  relatively 
unattenuated  through  the  medium.  In  this  range,  the  effective  medium  FE  models  give  a  very  good  approximation  to  the 
actual  medium  with  damage.  On  the  other  hand,  applying  the  homogenized  methods  in  the  range  with  input  wavelength 
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Fig.  16.  Output  signal  strength  v.s.  l//<a  for  tensile  waves  in  the  models  with  0o  =  0°.  The  crack  length  2 a  is  evenly  distributed  between  0.6  and  1.4,  with 
average  of  1.0.  Left  figure:  The  normalized  peak  averaged  stress  measured  at  distance  equal  to  10.  Right  figure:  The  normalized  peak  averaged  stress 
measured  at  distance  equal  to  20. 


below  the  “cut  off’  wavelength  is  inappropriate,  unless  we  modify  the  homogenized  methods  to  take  the  wave  length  depen¬ 
dent  energy  dissipation  into  account. 

We  also  expect  the  “cut  off’  1  / ka  to  decrease  when  the  crack  density  decreases.  This  is  because  when  the  crack  density 
decreases,  there  are  fewer  cracks  per  unit  area.  Therefore,  less  energy  will  be  dissipated  by  the  cracks  when  the  wave  passes 
through  the  media. 

On  a  PC  with  a  Pentium  4,  2.0  GHz  CPU  and  1GB  RAM,  it  took  from  3  to  30  h  for  a  single  point  in  the  tensile  Figs.  13-16 
depending  on  the  crack  distribution  and  the  time  frame  of  simulation  for  the  wave  propagation.  For  the  compressive  cases, 
the  time  varied  from  2  to  5  h. 

5.2.  Distributed  crack  lengths 

Fig.  16  shows  the  output  signal  strength  as  a  function  of  1  /ka  for  a  set  of  models  with  90  =  0°  and  average  crack  length  2a. 
These  are  different  than  the  previous  models  however,  because  the  individual  crack  length  is  variable.  For  simplicity,  we 
choose  an  even  distribution  of  crack  lengths  between  0.6  and  1.4  with  the  average  of  1.0,  (again,  except  for  the  cracks  that 
touch  the  boundaries).  Comparing  Fig.  13  with  Fig.  16,  we  can  see  that  for  the  models  with  distributed  crack  lengths,  the 
signal  strength  values  have  a  wider  band  for  the  range  of  10  <1  /ka<  60.  Outside  of  this  region,  the  results  in  Figs.  13  and 
16  are  quite  similar  to  each  other.  We  can  conclude  from  this  result,  that  a  distributed  crack  length  will  not  have  a  large 
effect  on  the  “cut  off’  wavelength  and  the  “safe  range”  for  using  the  homogenized  effective  medium  calculations. 

From  the  above  figures  we  can  see  that  in  most  of  the  cases  the  output  signal  strength  is  weaker  than  the  input  signal 
strength.  The  reason  is  obvious:  the  reflection  and  diffusion  between  the  cracks  reduce  the  average  signal  strength.  However, 
in  some  regions  of  several  samples,  we  see  that  the  output  signal  is  actually  stronger  than  the  input  signal.  This  can  be  ex¬ 
plained  by  the  random  distribution  of  the  cracks.  When  the  wave  travels  through  a  field  of  cracks  in  the  direct  numerical 
simulations,  part  of  the  wave  is  reflected  by  the  crack  surfaces  and  interferes  with  other  components  of  the  wave  which  fol¬ 
low.  Therefore,  the  local  distributions  of  the  cracks  in  the  models  can  actually  increase  the  local  stress  wave  amplitude  in 
some  regions. 

6.  Conclusions  and  discussion 

In  this  paper,  we  performed  direct  numerical  simulations  of  waves  traveling  in  micro-crack  damaged  media  and  com¬ 
pared  these  to  results  using  a  homogenized  effective  medium  calculation.  For  the  numerical  simulations,  we  built  several 
models  with  different  crack  distributions  for  each  of  three  angular  distributions;  90  =  0°,  90  =  30°  and  90  =  60°.  We  studied 
the  wave  speed  variation  as  a  function  of  the  input  wavelength  in  each  model  for  both  tensile  waves  and  compressive  waves. 
The  results  for  wave  speed  in  the  effective  medium  models  are  close  to  the  numerical  simulation  results. 

When  a  plane  wave  approaches  an  obstacle,  the  obstacle  will  reflect  (or  block)  a  portion  of  the  advancing  wave  energy.  If 
the  obstacle  is  small  with  respect  to  the  wavelength  of  the  advancing  wave,  there  will  be  a  small  disturbance  in  the  wave 
front,  but  the  majority  of  the  energy  will  pass.  Essentially,  the  wave  travels  around  the  obstacle.  However,  if  the  obstacle 
is  large  with  respect  to  the  wavelength,  it  will  reflect  a  significant  portion  of  the  wave,  detracting  from  the  energy  remaining 
in  the  transmitted  wave.  These  observations  are  supported  by  Eringen’s  [34]  analytical  solutions  of  scattering  of  a  plane 
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wave  with  a  single  cylindrical  obstacle.  Eringen  found  that  the  total  scattered  power  increases  with  decreasing  wavelength 
(i ka  =  2n a/X)  of  the  incident  plane  wave.  Here,  a  is  the  radius  of  the  cylinder  and  2  is  the  wavelength  of  the  incident  plane 
wave.  For  incident  waves  with  large  wavelengths  (large  1  / ka )  the  scattered  elastic  energy  approaches  zero  and  the  wave  does 
not  “see”  the  obstacle;  we  have  determined  that  for  crack  systems  modeled  using  our  finite  element  approach  this  occurs  for 
1  /ka  >  60  for  media  in  tension,  and  1  / ka  >  1 0  for  those  in  compression.  In  the  case  of  distributed  obstacles  (such  as  the  micro¬ 
cracks  in  the  present  study),  the  effect  is  more  evident  since  the  plane  wave  energy  that  passes  by  the  obstacles  (cracks),  will 
be  further  blocked  by  other  obstacles  behind  the  first  ones. 

We  also  studied  the  relationships  between  the  input  wavelength  and  output  signal  strength.  These  results  show  that  for 
both  tensile  and  compressive  waves,  there  is  a  “cut  off’  wavelength  above  which  the  wave  can  easily  pass  through  the  dam¬ 
aged  medium  and  that  this  “cut  off’  is  different  in  tension  and  compression.  When  the  input  wavelength  is  shorter  than  the 
“cut  off’  wave  length,  the  output  wave  strength  decreases  linearly  with  the  input  wavelength.  When  the  input  wavelength  is 
close  to  average  crack  length,  (l/ka  =  1)  the  wave  is  highly  diffused.  This  suggests  a  “safe  range”  for  using  the  homogenized 
effective  media  methods  in  studying  these  problems  without  considering  the  wave  length  dependent  energy  diffusion. 
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This  study  focuses  on  the  prediction  of  the  anisotropic  effective  elastic  moduli  of  a  solid 
containing  microcracks  with  an  arbitrary  degree  of  alignment  by  using  the  generalized 
self-consistent  method  (GSCM).  The  effective  elastic  moduli  pertaining  to  anti-plane  shear 
deformation  are  discussed  in  detail.  The  undamaged  solid  can  be  isotropic  as  well  as  aniso¬ 
tropic.  When  the  undamaged  solid  is  isotropic,  the  GSCM  can  be  realized  exactly.  When  the 
undamaged  solid  is  anisotropic  it  is  difficult  to  provide  an  analytical  solution  for  the  crack 
opening  displacement  to  be  used  in  the  GSCM,  thus  an  approximation  of  the  GSCM  is  pur¬ 
sued  in  this  case.  The  explicit  expressions  of  coupled  nonlinear  equations  for  the  unknown 
effective  moduli  are  obtained.  The  coupled  nonlinear  equations  are  easily  solved  through 
iteration. 

©  2009  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Microcracks  are  common  defects  in  solids  and  multiple  microcracks  usually  coexist  in  a  single  solid.  The  prediction  of  the 
effective  elastic  properties  of  a  microcracked  solid  is  technically  challenging  and  can  find  many  practical  applications. 
The  major  methods  developed  so  far  to  predict  the  effective  elastic  moduli  of  a  microcracked  solid  include  the  following: 
the  self-consistent  method  (SCM)  in  which  a  crack  is  embedded  directly  into  an  effective  medium  [1];  the  generalized 
self-consistent  method  (GSCM)  in  which  a  crack  is  surrounded  by  an  undamaged  matrix  region,  and  then  embedded  in 
the  effective  medium  [2,3];  the  Mori-Tanaka  method  (MTM)  [4];  the  differential  scheme  method  (DS)  [5-7];  and  the  mod¬ 
ified  differential  scheme  (MDS)  [8].  Recently  the  GSCM  in  conjunction  with  a  finite  element  method  (FEM)  was  developed  to 
take  into  account  crack  face  contact  and  friction  [9].  Most  recently  the  representative  unit  cell  approach  was  proposed  by 
Kushch  et  al.  [10]  to  calculate  effective  elastic  properties  of  a  microcracked  solid.  Here  it  shall  be  noted  that  the  effect  of  crack 
orientation  statistics  on  the  anisotropic  effective  moduli  of  a  microcracked  solid  was  in  fact  first  investigated  by  Santare  et  al. 
[3]  through  the  introduction  of  the  crack  orientation  distribution  function  </>($)  which  was  later  adopted  in  [11]  within  the 
framework  of  the  SCM  to  study  the  problem  of  cracks  with  an  arbitrarily  degree  of  alignment  in  a  material  that  is  originally 
anisotropic  before  the  damage  occurs. 

In  this  research  we  analytically  study  the  anisotropic  effective  elastic  moduli  of  a  solid  containing  microcracks  with  an 
arbitrary  degree  of  alignment  by  using  the  GSCM.  Our  model  is  in  principle  based  on  the  GSCM  developed  by  Santare 
et  al.  [3].  In  [3]  the  GSCM  was  used  to  predict  the  anisotropic  moduli  under  plane  stress  loading.  Here  we  are  concerned  with 
the  effective  elastic  moduli  pertaining  to  anti-plane  shear  deformation.  In  our  model  the  undamaged  material  can  be  isotro¬ 
pic  as  well  as  anisotropic.  An  exact  solution  to  the  cracked  elliptical  inclusion  problem,  which  is  essential  in  the  realization  of 
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GSCM,  can  be  derived  when  the  undamaged  material  is  isotropic.  On  the  other  hand  an  approximate  analytical  solution  can 
still  be  derived  when  the  undamaged  material  is  anisotropic  and  when  the  crack  density  is  not  very  high. 

2.  The  effective  moduli  of  a  microcracked  solid 


For  a  microcracked  solid,  the  strain  energy  relationship  between  the  effective  medium  and  the  actual  microcracked  solid 
can  be  expressed  as  [2,3] 


1 

2 


C*  /jO  (T®  — 

Jijklukluij  — 


lsmo°uo%+  1 


2V 


0) 


where  SJW  is  the  effective  compliance  of  the  damaged  material,  Sijki  is  the  compliance  of  the  undamaged  material  and  V  is  the 
sample  volume.  The  above  integral  is  the  energy  dissipated  through  the  opening  of  each  microcrack,  summed  over  all  M 
cracks,  and  erg  is  the  applied  homogeneous  stress  field  while  t?  =  is  the  traction  along  the  crack  face  if  the  crack  did 
not  exist  and  [u/]  is  the  crack  opening  displacement.  In  this  research  we  focus  on  the  two  dimensional  case  in  which  all 
the  cracks  penetrate  the  solid  through  the  x3-direction.  In  addition  we  only  discuss  the  effective  elastic  moduli  pertaining 
to  anti-plane  shear  deformation.  In  the  following  we  will  address  two  cases:  (i)  the  undamaged  material  is  isotropic;  (ii) 
the  undamaged  material  is  anisotropic. 


2.1.  Isotropic  undamaged  material 


The  degree  of  crack  alignment  can  be  described  by  the  crack  orientation  distribution  function  0(0)  with  0,  (|0|  <  n/2)  being 
the  angle  between  an  individual  crack  and  the  positive  Xi-axis  [3].  Without  losing  generality,  0(0)  can  be  taken  as  an  even 
function  of  9  since  we  have  assumed  that  the  undamaged  solid  is  isotropic.  For  simplicity,  0(0)  is  specifically  given  by  [3,9] 
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where  0O  <  n\2.  The  two  special  cases  of  perfectly  aligned  cracks  and  randomly  oriented  cracks  correspond  to  0O  =  0  and 
0O  =  7i/2,  respectively,  in  Eq.  (2). 

Once  we  have  introduced  0(0),  the  summation  in  the  energy  relationship  Eq.  (1)  can  be  written  as  an  integral  over  ori¬ 
entation  angle  0, 
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where  A  is  a  representative  area  of  the  sample. 

In  this  study,  we  assume  that  the  material  is  under  anti-plane  shear  deformation.  As  a  result  the  above  energy  relation¬ 
ship  Eq.  (3)  can  be  simplified  as 
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where  <7^  and  o\2  are  the  anti-plane,  far-held  stresses,  C44  and  C55  are  the  two  relevant  effective  moduli  of  the  damaged 
material,  p  is  the  shear  modulus  of  the  undamaged  material,  c  is  the  average  half  crack  length  and  p  =  Me2 \A  is  the  crack 
density  parameter. 

If  the  crack  did  not  exist,  the  uniform  traction  due  to  the  far-held  stresses,  t°,  along  the  line  of  the  crack  face  is  given  by 
t3  =  cos  00-33  -  sin  9g°3X  .  (5) 

Next  we  introduce  the  GSCM  [3]  to  approximately  take  into  account  the  interaction  between  the  cracks,  as  shown  in 
Fig.  1.  The  inclusion  is  assumed  to  have  the  properties  of  the  undamaged  material  with  known  elastic  moduli.  The  surround¬ 
ing  matrix  is  composed  of  the  effective  orthotropic,  damaged  material  with  the  principal  directions  along  the  x \  and  x2  axes. 
The  half-length  of  the  crack  is  c,  the  semi-major  and  semi-minor  axes  of  the  ellipse  are  a  and  b ,  respectively.  The  crack  den¬ 
sity  parameter  r\  relates  average  crack  length  to  the  area  of  the  ellipse,  but  in  general,  this  leaves  one  of  the  three  parameters 
a,  b  and  c,  unspecified.  Therefore,  as  an  additional  condition,  we  will  require  a2  -  b2  =  c2  to  be  satisfied.  This  is  the  same  rela¬ 
tionship  that  was  used  in  [3]  for  convenience,  but  here  it  is  necessary  in  order  to  make  analytical  solutions  possible.  The  two- 
phase  composite  is  subjected  to  uniform  anti-plane  shearing  o3A  and  o\2  at  infinity. 

By  using  the  complex  variable  method  [12,13],  the  crack  opening  displacement  for  the  elliptical  domain,  depicted  in 
Fig.  1,  [u3]  can  be  obtained  exactly  as 


[Us] 


4^c2  -  x2 

n\\  +  r-jr2(i  -r)] 


djU*  +  bC, 


(a  +  b)jw 


55  cos  Oat, 


ajr  +  bC44 

(a  +  b)ju * 


sin  Oat 


(1*1  <  c) 


(6) 
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where  a  =  ^(R  +  R  1 ),  b  =  § (K  -  ft  1 ),  (ft  >  1 ),  r  =  fi*lju  with  /T  =  ^/C44C55.  Here  the  parameter  R ,  which  is  introduced  during 
the  analysis  by  using  the  complex  variable  method  and  conformal  mapping,  can  be  expressed  in  terms  of  a  and  b  as  R  = 

A  detailed  derivation  can  be  found  in  the  Appendix. 

Consequently  the  integral  in  Eq.  (4)  can  also  be  carried  out  exactly  to  arrive  at  the  following  expression 


ay*  +  bQ4  7t[2e0  -  sin(2^0)](crg1)2  +  bC^  7t[20o  +  sin(2g0)]((Tjj2)2 

(a  +  b)n*  0O[1  +r-R-2(\ -T)]  +  (a  +  b)n*  0O[1  +r-R“2(l -  r)]\’ 

(7) 

which  immediately  leads  to  the  following  two  coupled  nonlinear  equations  for  C44  and  C55  by  observing  the  fact  that  and 
o%2  are  arbitrary 

^  |  a/i*  +  bC*44  nrj[290  -  sin(2 0o)] 

c;5“  +  (a  +  b)p  200 [1  +r-R-2(l  -T)}’ 

n  _  1  an"  +  bC)5  TU][29o  +  sin(2g0)] 

CU-  +  (a  +  b)n*  2e0[i  +r-r2(i  -r)]’ 

where  rj  =  ^  =  n{R2^R-2y  The  above  set  of  nonlinear  equations  can  be  solved  easily  through  iteration. 

In  the  following  we  discuss  in  more  detail  some  special  cases  for  the  above  solution. 


(^3°i)2  |  (^302)2  =  (^i)2  +  (^02)2  1  n 

2  C*55  ^  2  C44  2 fi  ijl 


2AA.  Aligned  cracks  (0o  =  0) 

In  the  case  of  aligned  cracks,  it  follows  from  Eq.  (8)  that 

C55  =  fa 

1  ar  +  b  2m]  (9) 

7*-  +  (a  +  bjr  i  +  r-R~2(i  -r)’ 

which  means  that  there  is  no  reduction  in  stiffness  along  the  Xi -direction  for  aligned  cracks. 

When  the  crack  density  parameter  is  extremely  low,  i.e.,  r]  <c  1,  Eq.  (9)2  further  reduces  to 

Q4  =  jU/(l  +  7T?7),  (10) 

which  is  exactly  the  non-interaction  approximation  (NIA)  for  aligned  cracks  described  by  Kachanov  [14]. 

Fig.  2  illustrates  C44  calculated  by  using  Eq.  (9).  The  dashed  line  is  the  result  of  NIA.  We  see  that  the  predictions  of  the 
GSCM  are  lower  than  that  of  NIA,  especially  when  r]  is  large. 


2 A .2.  Randomly  oriented  cracks  (60  =  n/2) 

In  the  case  of  randomly  oriented  cracks,  it  follows  from  Eq.  (8)  that 


X.  Wang  et  al.  / Engineering  Fracture  Mechanics  76  (2009)  1910-1919 


1913 


Fig.  2.  The  calculated  C*44  for  an  isotropic  solid  with  perfectly  aligned  cracks. 


through  which  r  can  be  exactly  determined  as 


r 


-( nrj  -2R  2)  +  \Jn2r]2  +4(1  -  nrjR  2) 

2(1  +R-2) 


(12) 


Eq.  (1 1  )i  shows  that  the  effective  material  properties  are  still  isotropic  when  the  cracks  are  randomly  oriented.  In  addition 
when  the  crack  density  parameter  is  extremely  low,  i.e.,  rj  <  1,  Eq.  (11)2  further  reduces  to 

ft  =  n/(  1  +  0.5711/),  (13) 

which  is  again,  the  NIA  described  by  Kachanov,  this  time  for  randomly  oriented  cracks  [14]. 

Fig.  3  shows  p  calculated  by  using  Eq.  ( 1 2 ).  The  dashed  line  is  the  result  of  NIA.  We  observe  a  similar  behavior  here  to  that 
of  C44  for  aligned  cracks.  We  also  observe  that  the  values  of  p  for  a  solid  with  randomly  oriented  cracks  are  higher  than  those 
of  C44  for  a  solid  with  aligned  cracks. 


2 A. 3.  An  arbitrary  degree  of  alignment  (0<  60<  n/2) 

The  present  model  can  be  used  conveniently  to  predict  the  anisotropic  effective  moduli  of  a  solid  containing  microcracks 
with  an  arbitrary  degree  of  alignment  through  the  introduction  of  the  crack  orientation  distribution  function  0(0). 

We  demonstrate  in  Figs.  4  and  5  the  predicted  values  of  C44  and  C55  as  functions  of  17  for  five  different  values  of  90  =  7i/16, 
71/6, 7i/4, 7i/3,  7t/2.  It  is  clear  from  the  two  figures  that  a  decrease  in  0O  will  result  in  a  decrease  in  C44  and  an  increase  in  C*55. 


Fig.  3.  The  calculated  ju*  =  C*44  =  Cg5  for  an  isotropic  solid  with  randomly  oriented  cracks. 
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Fig.  4.  The  predicted  values  of  C*44  for  a  cracked  isotropic  solid  as  a  function  of  r\  for  five  different  values  of  00- 


Fig.  5.  The  predicted  values  of  Cg5  for  a  cracked  isotropic  solid  as  a  function  of  r\  for  five  different  values  of  0O. 


The  anisotropy  of  the  effective  medium  monotonically  increases  as  0O  decreases  from  n\2  (for  randomly  oriented  cracks)  to 
zero  (for  perfectly  aligned  cracks). 

Due  to  the  fact  that  the  parameter  /i*  =  \/C*44C55  is  an  invariant  under  coordinate  rotation,  it  can  be  considered  as  a  mea¬ 
sure  of  the  overall  stiffness  of  the  cracked  material.  Therefore,  it  is  of  interest  to  check  how  p  is  influenced  by  0O  and  r\.  We 
demonstrate  in  Fig.  6  the  predicted  values  of  ju  as  a  function  of  p  for  four  different  values  of  0O.  It  is  observed  that  an  increase 
in  Oq  will  cause  a  decrease  in  the  overall  stiffness  characterized  by  fi\  especially  when  the  crack  density  rj  is  large.  This  fact 
can  be  more  clearly  observed  in  Fig.  7  in  which  ju  is  plotted  as  a  function  of  0O  for  six  different  values  of  r\.  Figs.  6  and  7  dem¬ 
onstrate  that  for  a  microcracked  solid,  randomly  oriented  cracks  (0O  =  7c/2)  gives  the  lowest  stiffness  among  all  possible  crack 
orientation  distributions,  while  aligned  cracks  (0O  =  0)  gives  the  greatest  stiffness. 

2.2.  Anisotropic  undamaged  material 


In  the  previous  discussion  we  have  assumed  that  the  undamaged  material  is  isotropic.  Next  we  consider  the  more  com¬ 
plex  situation  in  which  the  undamaged  material  is  anisotropic.  Without  losing  generality  we  can  assume  that  the  intact 
material  is  orthotropic  with  its  principal  directions  along  the  X\  and  x2  axes,  respectively.  When  the  intact  material  is  ortho¬ 
tropic,  we  cannot  necessarily  take  the  crack  orientation  distribution  function  </>(0)  as  an  even  function  of  6.  Here  </>(0)  takes 
the  following  form 


m 


02-01  5  ^1  <  0  —  ^2 
0,  6  >  02  or  6  <  0! 


(14) 
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Fig.  6.  The  predicted  values  of  /u*  for  a  cracked  isotropic  solid  as  a  function  of  rj  for  four  different  values  of  0o. 


Fig.  7.  The  predicted  values  of  //  for  a  cracked  isotropic  solid  as  a  function  of  60  for  six  different  values  of  rj. 


where  |0a|,  \02\  <  n  and  02  -  0\  <  n.  The  definition  region  of  6  in  Eq.  (14)  has  been  enlarged  relative  to  Eq.  (2),  to  \6\  <  n  to 
incorporate  the  complex  situation  in  which  the  average  orientation  of  the  cracks  9  =  (0i  +  e2)l 2  may  be  arbitrary  with  respect 
to  the  Xi,  x2  axes.  In  this  case,  perfectly  aligned  cracks  and  randomly  oriented  cracks  correspond  to  Q\  =  02  and  02  =  -0\  =  n\2 
(or  more  generally  but  equivalently  02-0i  =  n),  respectively,  in  Eq.  (14). 

Now  we  can  write  the  energy  relationship  in  Eq.  (3)  as  follows 


q4(cr°i)2  c;5oyg2  c;5(cr°2)2 

2  n*2  jr2  2 /i'2 


(<ii)2 ,  «)2 ,  n 

2  C55  2  C44  +  2c2 


f2m  j  [u3]t§, 

JB-i  JCk 


dCkd6 , 


(15) 


where  and  are  the  anti-plane,  far-held  stresses,  Q4,  and  C*45  are  the  three  relevant  effective  moduli  of  the  dam¬ 
aged  material,  /u*  =  \J Q4C55  -  C*45  >  0  and  C44  and  C55  are  the  two  material  moduli  of  the  undamaged  material.  When  the 

undamaged  material  is  anisotropic,  it  may  very  well  be  impossible  to  calculate  exact  values  for  [u3]  using  the  boundary-value 
problem  shown  in  Fig.  1.  Therefore,  in  the  following  we  use  an  approximation  to  the  BVP  to  calculate  [u3].  When  the  crack 
density  parameter  is  relatively  low,  we  have  R  >  1  or  equivalently  a  «  b  >  c.  As  a  result,  [u3]  can  be  obtained  approximately 
through  the  following  method:  (i)  first  calculate  the  uniform  stress  held  within  an  uncracked  circular  inclusion  with  material 
moduli  C44  and  C55  surrounded  by  the  effective  medium  with  material  moduli  C44,  C*55  and  C45  [15,16];  (ii)  then  solve  the 
problem  of  a  crack  in  an  inhnite  homogeneous  material  with  material  moduli  C44  and  C55  subjected  to  the  remote  loading 
which  is  equal  to  the  internal  uniform  stress  held  obtained  in  (i).  By  using  this  method  we  obtain  the  following  approximate 
expression  of  [u3] 
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N]  = 


2  Vc2  -  x2 


2  C, 


U** 


^55)a32  ^45°3\ 


[i  +  c 44  p*[i  +  r  +  p(i  —  r )] 


cos  a-  2'‘ 

/J.  +  C44  jU*[l  +  —  p(l  —  i^)] 


(M  <  C) 


(16) 


where  p  =  ^  /(  =  VQ4C55  and  r  =  fi/fi. 

y  c44+y  c55 

The  integral  in  Eq.  (15)  can  now  be  carried  out  to  arrive  at  the  following  relationship 

T0  /tO  f'*  f  /t-0  ^ 


Q4(^l)2 


2/r: 


^45^31  °32  |  ^55  (^32 ) 

. n  i 


2  fi*: 


2  C55 


+  ; 


nt] 


C44 

(/a*  +  c55)(g^2)  q5^30i^3°2 

2(Ju  +  C44) 

p*[i+r  +  p(i-r)] 

H  (Afl  +  ^44)  (^31 )  ^45*^31  ^32 

2(ju  +  c44)  /i*[i  +  r-p(i  -r)] 


«)2 

2C44  2/a(02  —  0i) 

[2(02  -  Oi)  +  sin(202)  -  sin(20!)] 

[2(02  -  0i )  -  sin(202)  +  sin(20i )] 


IA 

(fi*  +  q4)^3i^T32  qs((732)  ,  ^44 

(q  +  q5)^i^2  ^45(^31) 

[sin2  0t _ sin2  0il 

/A  +  C44 

jA*  [1  +  T”  —  p(l  —  T1)]  /A  +  C44 

p*[i+r  +  p(i-r)] 

3111  i/|  j 

(17) 

In  view  of  the  fact  that  a21  and  o\2  are  arbitrary,  we  arrive  at  the  following  three  coupled,  nonlinear  equations  for  Q4,  C*55 


and  C, 


44  _  j  ■  ^ 

'*2  C55  /i*(02  —  0l)(jU  +  C44) 


-55  . 
r  1*2 


-45 

,#*2 


1 


C44  /a *  (02  —  0i )  (p  +  C44) 


711/ 


(p*  +  q4)[2(02  -  0i )  -  sin(202)  +  sin(20i)]  |  C44q5[sin2  02  -  sin2  0i] 


2[\+r-p(\-r)\  p[i+r  +  p(i-r)]  _ 

C44(p*  +  q5)[2(02  -  0i )  +  sin(202)  -  sin(20i)]  C45[sin2  02  -  sin2  0i] 


2p[i+r  +  p(i-r)] 


i+r-p(i-r) 


(18) 


mj 


2p*(02  -0i)(p  +  C44)  _ 

fi*  +  C44 


C44C45 [2(02  -  01 )  +  sin(202)  -  sin(20i)]  |  C45[2(02  -  0Q  -  sin(202)  +  sin(20i)] 


+ 


2p[i+r+p(i-r)] 
C 44(11*  +  C55) 


2[i+r-p(i-r)] 


[sin2  02  -  sin2  0i] 


_i+r-p(i-r)  p[i+r+p(i-r)]_ 

which  can  also  be  solved  easily  through  iteration.  We  see  from  Eq.  (18)  that  when  </>(0),  defined  by  Eq.  (14),  is  an  even  func¬ 
tion  of  0,  i.e.,  0i  =  -02,  then  the  effective  medium  will  be  orthotropic  with  C45  =  0.  In  other  words,  when  the  average  orien¬ 
tation  of  the  cracks  is  zero  the  undamaged  material  and  the  effective  medium  will  possess  the  same  principal  axes.  In  the 
following  we  concentrate  our  discussion  on  the  effective  elastic  moduli  of  a  solid  containing  perfectly  aligned  cracks  and 
one  containing  randomly  oriented  cracks,  respectively. 


2.2 A.  Aligned  cracks  (6-1  =  02) 

When  the  cracks  are  perfectly  aligned,  we  have  0i  =  02  and  both  the  numerator  and  denominator  in  the  right  hand  side  of 
Eq.  (18)  vanish.  In  this  case  applying  the  L’Hospital’s  rule  to  Eq.  (18)  when  02  0i  yields  the  following: 


Q4=  1  nri 

p*2  C55  p*(p  +  C44) 


C44 


nrj 


(p*  +  C44)[l  -  cos (20!)]  C44C45  sin(20i) 


+  - 


i+r-p(i-r)  p[i+r  +  p(i-r)] 

C44(p*  +  C;5)[l  +  cos(20i)]  C45  sin(20i) 


p*(p  +  C44)  _  p[  1  +  r  +  p(  1  —  r )] 


i+r-p(i-r)J’ 


nrj 


C44q5[l +cos(20i)]  q5[l -cos(20i)] 


2p*(p  +  c44)  L p[i+r  +  p(i-r)]  i  +  r-p(i-r)  L i+r-p(i-r)  p[i  +r  +  p(i -r)]j 


/A* 


c44(p*  +  q5) 


sin(20i) 


(19) 

Now,  we  can  further  look  into  two  special  cases  of  the  above  formulation:  the  cracks  are  aligned  horizontally  (0i  =  02  =  0) 
and  the  cracks  are  aligned  vertically  (0i  =  02  =  7i/2).  When  0i  =  02  =  0,  it  follows  from  Eq.  (19)  that 


q5  —  0, 


C55  =  C55, 


1 


27i7/C44(/r 


c44  /a* ii(fi  +  C44)[1  +  r  +  p(i  —  r)]  5 

which  means  that  there  is  no  reduction  in  stiffness  along  the  Xi -direction  for  horizontally  aligned  cracks. 
On  the  other  hand  when  0i  =  02  =  7i/2,  it  follows  from  Eq.  (19)  that 


(20) 


Qs  =  0 


q4  =  c44, 

1  ,  2nrj(/A* 

C55 


q5  c55  p*(p  +  c44)[i  +  r-p(i  -r)y 


(21) 


which  means  that  there  is  no  reduction  in  stiffness  along  the  x2-direction  for  vertically  aligned  cracks. 
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Fig.  8.  The  variations  of  C44  (O-i  =  02  =  0)  and  C55  (0i  =  02  =  n/2)  for  an  orthotropic  solid  C44  - 1  .2 n  and  C55  =  0.8333 /i  containing  perfectly  aligned  cracks  as 
functions  of  rj. 


Fig.  9.  The  variations  of  C*44,  C*55  and  C45  as  functions  of  6-1  (=02)  for  an  orthotropic  solid  C44  =1.2/1  and  C55  =  0.8333 /i  containing  perfectly  aligned  cracks  with 
rj=lln. 


Fig.  8  illustrates  the  variations  of  C44  for  0i  =  02  =  0  and  C55  for  0a  =  02  =  n/2  as  functions  of  y\.  To  plot  the  results  on  one 
graph,  the  undamaged  anisotropic  moduli  are  specifically  chosen  as  C44  =  1.2/7  and  C55  =  0.83  33 /7.  It  is  interesting  to  observe 
that  C44  for  =  02  =  0  and  C*55  for  6\  =  02  =  n/2,  both  of  which  are  decreasing  functions  of  77,  converge  to  basically  the  same 
value  as  r\  1.  We  also  observe  that  the  effective  material  is  isotropic  (C44  =  C*55  =  0.8333/7)  at  low  crack  density  r\  «  0.10  for 
horizontally  aligned  cracks  (0i  =  02  =  0). 

We  can  also  see  from  Eq.  (19)  that  when  the  cracks  are  not  aligned  horizontally  or  vertically,  the  principal  directions  of 
the  anisotropic  effective  medium  are  no  longer  the  same  as  those  of  the  undamaged  material.  We  illustrate  in  Fig.  9,  the  vari¬ 
ations  of  C44,  C^5  and  C*45  as  functions  of  0a  (=02)  for  an  orthotropic  solid  with  undamaged  moduli,  C44  =  1.2^i  and 
C55  =  0.8333^  containing  aligned  cracks  with  ^  =  1/71.  We  observe  from  Fig.  9  that:  (i)  C45  attains  a  maximum  value  of 
C*5  =  0.3073//  at  0^  =48°  (the  reason  why  C*45  does  not  attain  its  maximum  value  at  45°  is  due  to  the  anisotropy  of  the 
undamaged  material);  (ii)  C44  is  an  increasing  function  and  C55  is  a  decreasing  function  of  0i.  As  expected 
C44  =  C44  =  1.2/7  for  vertically  aligned  cracks  (0i  =  90°)  and  C*55  =  C55  =  0.8/7  for  horizontally  aligned  cracks  (0t  =  0).  More 
interestingly,  our  numerical  results  also  verify  that  the  effective  compliance  constants  Sj5  =  -C*45/fi*2  and  S*55  =  Q4//7*2 
in  the  Cartesian  coordinate  system  (x\,x'2),  which  is  rotated  by  the  crack  alignment  angle  0i  with  respect  to  the  principal 
coordinate  system  (xi,  x2),  are  exactly  the  same  as  the  compliance  constants  S45  =  —  C45// 72  and  5 55  =  C44//72  for  the  undam¬ 
aged  material  in  the  same  coordinate  system  (x\,x'2)  . 
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Fig.  10.  The  variations  of  C*44  and  C55  for  an  orthotropic  solid  C44  =  1.2/r  and  C55  =  0.8333 /i  containing  randomly  oriented  cracks  as  functions  of  r\. 


2.2.2.  Randomly  oriented  cracks  (02  =  -0j  =  n/2) 

When  the  cracks  are  randomly  oriented,  Eq.  (18)  reduces  to 

C45  -  0, 

1  =  1  1  nrijfi'  +  Cy) 

c;5  C55  li^ii  +  c44)  [i  +  r  -  p(i  -  r)} 5  (22) 

1  =  1  nr}CAA(ii* +  C55) 

C44  C44  +  C44)  [i  +  /"”  +  p(i  —  r )] 

From  the  above  expressions,  we  can  see  that  C44  ^  C*55  when  C44  7^  C55  even  when  the  cracks  are  randomly  oriented.  To 
illustrate  this  fact  more  clearly  we  present  in  Fig.  10,  the  variations  of  C44  and  C*55  for  an  orthotropic  solid  containing  ran¬ 
domly  oriented  cracks  as  functions  of  rj.  The  undamaged  anisotropic  moduli  are  chosen  as  C44  =  1.2 p  and  C55  =  0.83  33 p.  It 
is  clear  that  C44  >c;5,  and  in  addition  the  difference  in  C44  and  Q5  decreases  as  the  crack  density  increases. 

Finally  it  is  of  interest  to  point  out  that  when  the  undamaged  material  is  isotropic ,  our  calculations  show  that  the  result  by 
using  the  approximation  of  the  GSCM  always  lies  between  that  by  using  the  exact  GSCM  and  that  by  using  NIA.  In  view  of  the 
fact  that  the  GSCM  is  only  an  approximate  scheme  to  account  for  crack  interactions,  the  predictions  derived  by  using  the 
GSCM  are  not  necessarily  more  accurate  than  those  derived  by  using  an  approximation  to  the  GSCM. 


3.  Conclusions 

By  means  of  the  GSCM,  we  analytically  investigated  the  anisotropic  effective  moduli  of  a  cracked  solid  subjected  to  anti¬ 
plane  shear  deformation.  When  the  undamaged  solid  is  isotropic,  the  two  coupled,  nonlinear  equations  for  the  two  effective 
moduli  C44  and  C*55  were  obtained  and  shown  in  Eq.  (8).  When  the  undamaged  solid  is  anisotropic,  the  three  coupled  non¬ 
linear  equations  for  the  three  effective  moduli  Q4*  C55  and  C45  were  derived  and  shown  in  Eq.  (18).  Detailed  numerical  re¬ 
sults  were  presented  to  illustrate  how  the  anisotropic  effective  moduli  are  influenced  by  various  degrees  of  crack  alignment 
(characterized  by  0O  for  the  isotropic  undamaged  material  or  0i  and  02  for  the  anisotropic  undamaged  material)  and  crack 
density  (characterized  by  rj). 
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Appendix  A.  Detailed  derivation  of  Eq.  (6) 

It  is  convenient  to  solve  the  boundary-value  problem  in  the  Cartesian  coordinate  system  (x,y,  z)  with  z  =  x3,  as  shown  in 
Fig.  1 .  First  the  remote  uniform  loading  o°zx  and  azy  can  be  written  as 
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G'L  =  (To-.  COS  t 


(Jo9  sin  t 


y  2  ^  '-vj  1/  i  iy  22 

(7°y  =  0-32  cos  0  -  (7^  sin  t 
Secondly  the  effective  moduli  of  the  damaged  material  in  the  (x,  y,  z)  coordinate  system,  Q, 


(Al) 


44  ’ 


C*55  and  CJ5  can  be  derived 


as 


C55  =  C55  cos2  6  +  Q4  sin2  0, 


C44  =  C44  cos2  < 


-  c;5  sin 


(A2) 


C45  =  (C44  -  Cg5)  cos  0  sin  { 


Due  to  the  fact  that  the  cracked  elliptical  inclusion  is  isotropic,  if  one  is  only  concerned  with  the  field  within  the  inclusion, 
it  is  sufficient  to  treat  the  matrix  as  isotropic  with  shear  modulus  p  subject  to  virtual  remote  uniform  shear  stresses  <7°x  and 
d°zy  such  that  [13] 

(a  -  ipb)(a°zx  +  pa°zy) 


<  ~  i<  =  - 


(a  +  b)Im{p} 
where  the  complex  constant  p  is  given  by 


P  = 


C45  +  ip*  _  (C55  C44)  cos  9 sin  9  +  i  \/C44 C5: 


C44  cos2 1 


■  Qs  sin 


(A3) 


(A4) 


Here  we  are  only  interested  in  the  expression  for  a°y  since  the  loading  a°x  will  not  induce  any  crack  opening  displacement. 
It  follows  from  Eqs.  (Al),  (A3),  and  (A4)  that 

■|2l"°  bRe{p}<7°2  -  [Qlm{p}  +  b\pf]a^ 


do  bRe{p}o~°i  +  [Qlm{p}  +  b|ph<rjj2  cps  g 


(a  +  b)Im{p} 


(a  +  b)Im{p} 


sind 


a/r  +  bC5 , 

(a  +  b)p* 


cos  9G32  ~ 


a/T  +  bC4, 

(a  +  b)p- 


sin  Pa?, . 


(A5) 


On  the  other  hand  it  can  be  deduced  from  Ref.  [12]  that  the  crack  opening  displacement  [u3]  due  to  remote  uniform  load¬ 
ing  <r°y  can  be  expressed  as 


M  =  ■ 


Aa^Vc2  -x2 


(\X\  <  c). 


p{\  +r-R~2(\  -r)\ 

Starting  with  Eqs.  (A5)  and  (A6),  we  can  easily  arrive  at  Eq.  (6). 


(A6) 
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Abstract  In  this  study  we  first  obtain  the  explicit 
expressions  for  the  15  effective  reduced  elastic  compli¬ 
ances  of  an  elastically  anisotropic  solid  containing  mul¬ 
tiple  microcracks  with  an  arbitrary  degree  of  alignment 
under  two-di  mensi  onal  def ormati  ons  w  i thi  n  the  frame¬ 
work  of  the  non-interaction  approximation  (NIA). 
U  nder  special  situations,  our  results  can  reduce  to  the 
classical  ones  derived  by  Bristow  (J  Appl  Phys  11: 
81-85, 1960),  and  M  augeand  Kachanov  (J  M  ech  Phys 
Solids  42(4):561-584,  1994).  Some  interesting  phe¬ 
nomena  are  also  observed.  For  example,  when  the 
undamaged  solid  is  orthotropic,  the  effective  in-plane 
shear  modulus  is  dependent  on  the  degree  of  the  crack 
alignment.  The  NIA  method  is  then  extended  to  obtain 
the  effective  electroelastic  properties  of  an  anisotropic 
piezoelectric  solid  containing  two-dimensional  insulat¬ 
ing,  permeable  or  conducting  microcracks  with  an 
arbitrary  degree  of  alignment.  We  also  derive  a  set  of 
fifteen  coupled  nonlinear  equations  for  the  unknown 
effective  reduced  elastic  compliances  of  a  microcrac¬ 
ked,  anisotropic,  elastic  solid  by  using  the  generalized 
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self-consistent  method  (GSCM).  The  set  of  coupled 
nonlinear  equations  can  be  solved  through  iteration. 

Keywords  Non-interaction  approximation  • 
Effective  properties  •  M  icrocracks  •  Elastic  material  • 
Piezoelectric  material 


1  Introduction 

The  prediction  of  the  effective  elastic  properties  of 
anisotropic  materials  containing  microcracks  has  been 
a  topic  of  micromechanics  for  nearly  30years  (see  for 
example,  Gottesman  etal.  1980;  Hashin  1988;  Mauge 
and  K  achanov  1994;  F  eltman  and  Santare  1999;  among 
others).  The  limitations  of  the  previous  works  lie  in  the 
following:  (i)  the  two-dimensional,  anisotropic  matrix 
studied  so  far  is  confined  to  the  special  case  in  which 
the  in-plane  deformations  and  the  out-of-plane  defor¬ 
mations  are  decoupled  (see  for  example,  M  auge  and 
Kachanov  1994);  (ii)  in  most  of  the  studies  the  micro¬ 
cracks  were  assumed  to  be  either  perfectly  aligned  or 
randomly  oriented  (M  auge  and  Kachanov  1994). 

In  this  study  we  first  address  the  problem  of  the 
effective  elastic  properties  of  a  generally  anisotropic 
elastic  matrix  containing  microcracks  with  an  arbitrary 
degree  of  alignment  under  two-dimensional  deforma¬ 
tions.  Here  we  consider  any  degree  of  crack  alignment, 
from  perfectly  aligned  to  completely  random,  through 
the  introduction  of  a  crack  orientation  distribution  func¬ 
tion  (Santare  etal.  1995).  Very  concise  expressions  for 
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the  effective  reduced  el  asti  c  compl  i  ances  of  the  cracked 
material  are  obtained  by  making  use  of  energy  con¬ 
siderations  and  by  ignoring  crack  interactions.  During 
the  development,  we  adopt  some  results  (for  example, 
the  crack  opening  displacement)  based  on  the  Stroh 
formalism  for  two-dimensional  anisotropic  elasticity 
(Ting  1996).  The  current  result  based  on  a  non-interac¬ 
tion  approximation  (NIA)  still  retains  accuracy  at  sub¬ 
stantially  higher  crack  densities  due  to  the  cancella¬ 
tion  of  the  competing  interaction  effects  of  shielding 
and  antishielding  (Kachanov  1992;  M  auge and  Kacha¬ 
nov  1994;  Kachanov  2007).  The  NIA  method  is  then 
extended  to  study  the  effective  electroelastic  proper¬ 
ties  of  a  microcracked,  anisotropic,  piezoelectric  solid. 
Here  the  microcracks  can  be  insulating  (Pak  1990;  Suo 
etal.  1992),  permeable(Suo  etal.  1992)  or  conducting 
(M  cM  eeking  1987;  Suo  1993).  It  is  of  interest  to  point 
out  that  the  non-convex  electric  enthalpy  used  for  insu¬ 
lating  cracks  is  different  than  the  convex  energy  den¬ 
sity  function  used  for  conducting  cracks.  The  effective 
elastic  moduli  of  a  microcracked,  generally  anisotropic 
elastic  solid  are  also  discussed  within  the  framework  of 
th  e  g  en  eral  i  z ed  sel  f -  c o  n  si  sten t  m  eth o d  ( G  S  C  M  ) ,  which 
approximately  takes  into  account  crack  interaction 
(Aboudi  and  Benveniste  1987;  Santare  et  al.  1995). 
It  is  found  that  some  observations  made  by  using  the 
NIA  assumption  are  no  longer  valid  when  using  the 
GSCM . 


2  The  theory  within  the  framework  of  NIA 

T  he  strai  n  energy  bal ance  equati  on  for  a  mi crocracked, 
elastically  anisotropic  solid  can  be  expressed  as 
(E shelby  1956;  Benveniste  1985;  Aboudi  and  Benven¬ 
iste  1987) 

l  l  1  M  r 

2Sijklaklaij=2SijkiaklaU+2^  X  / 

k=lCk 

(1) 

where  S*jkl  is  the  effective  compliance  of  the  dam¬ 
aged  material,  is  the  compliance  of  the  undamaged 
material,  V  is  the  sample  volume,  M  is  the  number  of 
microcracks  within  the  solid,  erP.  is  the  applied  homo- 
geneous  stress  field  while  tf  =  a^n,  is  the  traction 
along  the  crack  face  if  the  crack  did  not  exist  and  [«,-] 
is  the  crack  opening  displacement. 


Fig.  1  A  crack  of  half  crack  length  a  oriented  at  an  angle  9  to 
the  positive  xi-axis 


In  this  study  we  assume  that  all  the  cracks  pene¬ 
trate  the  solid  through  thex3-axis.  The  degree  of  crack 
alignment  can  be  described  by  the  crack  orientation 
distribution  function  <A(0)  with  0,  (|0|  <  n/2)  being 
the  angle  between  an  individual  crack  and  the  positive 
xi-axisas shown  in  Fig.  1  (Santare etal.  1995).  Without 
losing  generality,  4>(9)  can  be  taken  as  an  even  func¬ 
tion  of  0  such  that  the  average  orientation  is  parallel 
tothexi-axis.  For  simplicity,  4>(9)  is  specifically  given 
by  (Santare etal.  1995;  Wang  etal.  2009) 


0(0)  = 


2^’ 

0,  |0|  >  0o 


(2) 


where  0o  <  n/ 2.  The  two  special  cases  of  perfectly 
aligned  cracks  and  randomly  oriented  cracks  corre¬ 
spond  to  0o  =  0  and  0o  =  n/2,  respectively  in  Eq.  2. 

Asa  result,  under  two-dimensional  deformations  in 
which  the  displacements  depend  on  xi  and  xi  only, 
Eq.  1  can  be  further  expressed  as  (Santare  et  al.  1995; 
Ting  1996) 

i(<r0)rS'V0  =  ^(<r°)7'S'<r0 

+Ia  j  H9) /  (tO)TMdc*d0’ 

-So  Ck 

(3) 


where  A  is  the  sample  area,  and 

_o  r  o  ^.o  _o  _o  o]r  /d\ 

a  -  ^crn  a22  er23  cr13  al2 J  ,  (4) 
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[u]  = 

[[Ml]  [M2]  [M3]]7  . 

(7) 

We  add  that  =  Sap-Sa3S3p/S33  in  Eq.  5  are  the 
reduced  elastic  compliances  of  the  undamaged  solid, 
and  sap  =  saft  -  sa3stp/sh  in  E9- 6  are  the  effective 
reduced  elastic  compliances  of  the  cracked  solid  (Ting 
1996). 


The  traction  vector  t°  is  given  by 


t°  = 


r'Tl°2 
r0 


cose 


sin  e 


ct22  cos  e  -  a12  si  n  e 
i_ct2°3  cose  —  o' 33  sine 


(8) 


If  we  ignore  the  crack  interactions,  i.e.,  ??  <1,  then  the 
crack  opening  displacement  vector  [u]  =  [[mi]  [m2] 
[m3]]t  can  be  simply  given  by  (Ting  1996) 


[u]  =  2 Va2  —  x2L  H0,  (|x|  <  a)  (9) 


where  a  is  the  half  crack  length,  the  3  x  3  real  and 
symmetric  B  arnett- L  othe  tensor  L  for  the  undamaged 
solid  is  positive  definite  and  L_1  is  further  written  in 
the  foil  owing  form  for  the  convenience  of  thefol  lowing 
theoretical  development 


L  1 


Tn  Tn  Tib 
T12  T22  T23 

Tib  t23  y33 


(10) 


We  stress  here  that  in  writing  Eq.  9  we  have  employed 
the  property  of  the  second-order  tensor  L  under  coor- 
dinate  rotation  (also  see  Fig.  114  on  page  421  in  Ting 
1996). 

Substituting  Eqs.2,  8  and  9  into  Eq.  3,  we  can  arrive 
at 


2«r  )  S  <,  =j(a  )  S, 


a 

00 

/  V —  x26x 

J 

[  (t°)rL-1t°d0 

J 

_-6»° 

where  rj  =  Ma2/A  is  the  crack  density  parameter. 

We  can  clearly  observe  from  the  above  expression 
that  S'*  is  also  positive  definite  due  to  the  fact  that 
the  value  of  the  second  term  on  the  right  hand  side  of 
Eq.  11  is  always  positive  (i.e.,  the  elastic  energy  density 
is  always  increased  through  the  introduction  of  the  mi¬ 
crocracks).  The  two  integrals  in  Eq.  11  can  be  exactly 
carried  out  as 

J  sfa 2  —  x2d.v  =  — — ,  (12) 

—a 

and 

0o 

J  (t°)rL"1t°d0  =  ^[20o  -  sin(20o)](oii)2 

-00 

+^[20o  +  sin(20o)](a2°2)2 

+^[20o  +  sin(20o)](CT2O3)2 

+^[20o-sin(20o)](CTiO3)2 

+^{Tn[20o  +  sin(20o)] 
+T22[20O  -  sin(20o)]}  (of2)2 
+Fi3[20o  -  sin(20o)]a1°1ai°3 
+Ti2[20o  -  sin(20o)]a1Vi°2 

+T23[20O  +  sin(20o)]CT2o2CT2°3 
+Ti2[20o  +  si  n  (2f?o)]cr22cr  12 
+Tib[20o  +  Si  n  (2£?o)]cr23cr12 
+T23[20O  -  sin(20o)]a1ViO2- 

(13) 

In  view  of  the  fact  that  the  applied  homogeneous 
stresses  er-^,  ct2°2,  cr2°3,  <xi°3  and  (t1°2  can  be  arbitrary, 
then  we  obtain  the  following  expressions  for  the  effec¬ 
tive  reduced  elastic  compliances 
*  TC  nY'\'\ 

5 ii  =  S'n  +  -  sin(26»0)], 

4  =  4  +  ^j^[20o  +  sin(20o)L 

si 4  =  s'4  4  +  ^^[20o  +  sin(20o)L 
4^o 

4  =  ^5  +  ^^[20o-sin(20o)], 

4  =  See  +  (Tii[20o  +  sin(20o)] 

+T22[20o  -  sin(20o)]) , 
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4 2  =  S'n, 

Su  =  4’ 

4s  =  4  +  ^>0  -  sin(26»o)], 

*  tctiYm 

4  =  4  +  -^[20o  -  sin(20b)], 

4  =  4  +  +  sin(20o)L 

4  =  4. 

4  =  4  +  +  sin(26?0)], 

4  =  4’ 

4  =  4  +  +  sin(26>o)], 

4  =  4  +  ^[20o  -  sin(20o)].  (14) 

Equation  14  indicates  that  the  four  compliance  compo¬ 
nents  4’  44’  4s  apd  4s are  eclual to  correspond¬ 
ing  ones  for  the  undamaged  material  and  that  4  > 
S’n,  (i  =  1,  2, 4,  5,  6)  since  Yu,  Yu,  T33  >  0.  It  is  of 
interest  to  notice  that  the  present  result  of  4  =  4  is 
consistent  with  the  observation  by  M  auge  and  Kacha¬ 
nov  (1994)  that  vf2/Ef  =  vn/Ei  for  an  orthotropic 
solid  containing  perfectly  aligned  or  randomly  oriented 
cracks  under  plane  stress  deformation.  The  result  we 
obtain  is  valid  for  generally  anisotropic  materials  con¬ 
taining  cracks  with  an  arbitrary  degree  of  alignment. 

When  the  undamaged  material  is  orthotropic,  it 
follows  from  the  above  expression  that 

1  1  7tr] 

H*l2  11 12  400 

'200  +  sin  (2 0O)  200  -  sin  (20o)' 

L  EU  L22  J 

(15) 


which  clearly  demonstrates  that  the  effective  in-plane 
shear  modulus  n\2  depends  on  the  degree  of  crack 
alignment  characterized  by  0o  when  L\\  ^  L22.  Fur¬ 
thermore  when  the  undamaged  material  is  isotropic 
such  that  L\\  =  L22  =  with  n  and  v  being 
the  shear  modulus  and  the  Poisson's  ratio,  respectively, 
then  Eq.  15  reduces  to 


* 

I1 12 


1  +7ir](l  -  v )  ’ 


(16) 


which  is  the  classical  NIA  result  obtained  by  Bristow 
( 1960),  and  i s  i ndependent  of  the  degree  of  crack  align¬ 
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Fig.  2  The  variations  of  n{2  of  a  cracked  orthotropic  elastic 
material  as  a  function  of  Oo  by  using  Eq.  15  for  three  different 
values  of  Ln/Lji  =  0.5, 1.0, 1.5  with  nruin/Ln  =  0.25 


ment  (Kachanov  1992).  We  illustrate  in  Fig.  2  the  var¬ 
iation  of  n\2  as  a  function  of  0o  by  using  Eq.  15  for 
three  different  values  of  L11/L22  =  0.5, 1.0, 1.5  with 
TTrjun/Ln  =  0.25.  Weobservefrom  Fig.  2  that /xf2  is 
an  increasing  function  of  0o  when  Ln  <  L22;  whereas 
it  is  a  decreasing  function  of  0o  when  Ln  >  L22. 

In  the  foil  owing  we  discuss  two  typical  cases  for  the 
crack  orientations: 


2.1  Perfectly  aligned  cracks 

When  the  microcracks  are  perfectly  aligned,  we  have 
0o  =  0.  Consequently  it  follows  from  Eq.14  that 

4  =  4. 

4  =  S22  +  7tr]Y22, 

*!>44  =  S44  +  nriYx, 

4  =  4’ 

4  =  4  +  ni)Yll, 

4z  =  4> 

4  =  4. 

*4  =  4’ 

4  =  4. 

4  =  4  +  7T??T23> 

“4  =  4> 

$26  =  *4  +  xr)Y\2, 
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S4b  -  ^45’ 

^46  =  ^46  +  ^  ^^13> 

4  =  4- 


(17) 


In  addition,  if  the  undamaged  material  is  orthotropic, 
E  q.  17  further  reduces  to  E  q.  27  derived  by  M  auge  and 
Kachanov  (1994)  in  view  of  the  fact  that  Yu  and  Y22 
for  orthotropic  materials  are  explicitly  given  by  (Suo 
1990) 


Yu 

Y22 


1 


(18) 


Here  it  should  be  noticed  that  M  auge  and  Kachanov 
(1994)  discussed  the  plane  stress  problem. 


2.2  Randomly  oriented  cracks 


In  addition  if  the  undamaged  material  is  orthotropic, 
Eq.  19  reduces  to  Eq.32  derived  by  M  auge  and  Kacha¬ 
nov  (1994).  It  seems  that  the  factor  'jt/4'  appearing  in 
the  third  expression  of  Eq.32  in  Mauge  and  Kachanov 
(1994)  should  read  ‘n/2’ .  Furthermore  it  follows  from 
E  q.  19  that  the  overal  I  effective  el  asti  c  properti  es  w  i  1 1  be 
transversely  isotropic  for  an  originally  isotropic  mate¬ 
rial  containing  randomly  oriented  cracks  (0o  =  jt/2). 

Once  the  effective  reduced  elastic  compliances  have 
been  obtained,  the  effective  stiffnesses  can  be  conve¬ 
niently  obtained  as 


C°*  =  (S'*)-1 

where 


rcfi 

4*2 

C14 

C15 

4*6 1 

4*2 

4* 2 

C24 

C25 

4*6 

c14 

l24 

c44 

l45 

l46 

4 

4*s 

c45 

4*s 

4*6 

1-4*6 

4*6 

c46 

4*6 

4*6  J 

(20) 


(21) 


W  hen  the  microcracks  are  oriented  randomly,  we  have 
do  =  7t/2.  Consequently  it  follows  from  Eq.  14  that 

4  =  5 11  + 

,*  ,  1 
S22  —  S22  +  2nrlY22> 

»  1 

4  =  S44  +  27TrlY'&’ 

4  =  4  +  2 71  ^ 5/33  ’ 

4  =  4  +  \*i(Yu  +  y22), 

4  =  4. 

4  =  4. 

4  =  4  +  27T??^13’ 

4  =  4  +  2nilYl2, 

»  1 

4  =  4  +  27TrlY23 , 

4  ~  4> 

4  =  4  +  27T0^12’ 

4  =  4> 

*  1 

4  =  4  +  27r^^13, 

4  =  4  +  \xnYn.  (19) 


3  Extension  to  piezoelectric  solids 

In  the  previous  section  we  have  derived  the  effective 
reduced  elastic  compliances  of  a  purely  elastic  solid 
containing  microcracks.  In  fact,  the  energy  method  can 
also  be  conveniently  adopted  to  derive  the  effective 
electroelastic  properties  of  a  piezoelectric  solid  con¬ 
taining  insulating,  permeable  or  conducting  micro¬ 
cracks. 


3.1  A  piezoelectric  solid  containing  insulating 
microcracks 

In  this  case,  the  energy  balance  relationship  for  two- 
dimensional  problems,  in  which  the  displacements  ut 
and  the  electric  potential  <f>  depend  on  x\  and  X2  only, 
can  be  expressed  i  n  the  f ol  I ow i  ng  f orm  w hen  the  pi ezo- 
electric  solid  containing  insulating  microcracks  is  sub¬ 
jected  to  the  homogeneous  stresses  <rP.  and  the  homo- 
geneous  electric  displacements  £>9  (Suo  etal.  1992) 

'40\7Vr'0*\-1^0  '*'/'^0\7Vr'0\-1^0 

-(<t  )  (C  )  a  =  -(a  )  (L  )  a 

e0 

-00  Ck 

(22) 
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where 

~o  r  o  o  o  o  o 
a  -  L^ll  or22  or23  a13  (7 12 

[U]  =  [[«!]  [u2\  [mb]  [<p]]T, 


c°  = 


D\ 


D°2]T , 


(23) 

(24) 


c°*  = 


C11 

c  12 

C14 

c  15 

Ci6 

e\i 

£21  " 

C\2 

C22 

C24 

C  25 

C26 

£12 

£22 

C14 

C24 

C44 

C45 

C46 

£14 

£24 

Cl5 

C25 

C45 

C55 

C56 

eib 

^25 

5 

Cl6 

C26 

C46 

c% 

C66 

^16 

^26 

en 

£12 

ei4 

ei5 

616 

-€ll 

-*12 

.  e2\ 

£22 

«24 

£25 

626 

-*12 

~*22_ 

( 

[2 

r^ii 

r<* 

C12 

c14 

l15 

C16 

ell 

e21 

~ 

ch 

r<* 

c22 

c24 

C’25 

C26 

el2 

e22 

c14 

r<* 

C24 

c44 

u45 

c46 

e\4 

e*24 

cTs 

C* 

C25 

c45 

C5*5 

C* 

C56 

e15 

e25 

^26 

C46 

C5*6 

^66 

e16 

e26 

en 

eh 

e14 

e15 

e16 

_6n 

~612 

-  e21 

e22 

e24 

e25 

e26 

-^1*2 

-£22 

_ 

rCT1°2  cos  e  -  orfi  sin  O' 


t°  = 


cr22  COS  0 

cr2°3  COS  6 


a{2  sin  9 
cr-^  si  n  9 


D§  COS0  -  D?  sin  0 


(26) 


(27) 


_  Ls 2  IUj  U  —  LJ ^ 

In  Eqs.25  and  26,  C,j .  el}  and  e,2  are  the  elastic,  pie¬ 
zoelectric  and  dielectric  constants  of  the  undamaged 
solid,  while  C*. ,  e*  and  e*  are  the  corresponding 
(unknow  n)  constants  for  the  damaged  sol  i  d.  T  he  energy 
expression  in  Eq.22  is  in  fact  the  electric  enthalpy 
(Suo  et  al.  1992).  When  ignoring  crack  interactions 
and  assuming  that  the  crack  surfaces  are  traction-free 
and  charge-free  (Pak  1990;  Suo  etal.  1992),  the  jump 
in  displacement  and  electric  potential  across  the  crack 
surface  can  be  simply  given  by 


[u]  =  2  \/  m2  —  x2L  1t° 


<  a) 


(28) 


where  L  is  the  4  x  4  real  and  symmetric  Barnett- 
Lothe  tensor  for  the  undamaged  piezoelectric  solid 
(Ting  1996).  We  add  that  L  is  not  positive  definite  and 
can  be  further  expressed  as 
Lll  L14 


L  = 


T  T 

.14 


-L44 


(29) 


where  the  3  x  3  symmetric  matrix  Ln  is  positive  def¬ 
inite  and  L44  >  0. 

Substituting  the  above  results  into  Eq.  22,  we  finally 
arrive  at  a  concise  expression  for  the  effective  elec¬ 
troelastic  properties  of  the  microcracked  piezoelectric 
solid 
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(c°*)-i  =  (cv1 + 

+  „,[20o  +  sin(20,)]  ,h  (30) 
40q 


where 


“1  0  0  0  0  0  0" 

0  0  0  0  1  0  0 

0  0  0  1  0  0  0 

_0  0  0  0  0  1  0_ 

“0  0  0  0  1  0  0“ 

0  1  0  0  0  0  0 

0  0  1  0  0  0  0 

0  0  0  0  0  0  1 


(31) 


It  can  be  easily  verified  from  Eqs.  29  and  30  that  both 

are  positive  def¬ 


C°*  defined  by  Eq.  21  and 


:11  fc12 


inite.  It  can  also  be  observed  from  Eq.  30  that  the  eigh¬ 
teen  components  (12)  =  (21),  (13)  =  (31),  (24)  =  (42), 
(34)  =  (43),  (26)  =  (62),  (36)  =  (63),  (17)  =  (71),  (47)  = 
(74)  and  (67)  =  (76)  of  (C0*)-1  do  not  change  due  to 
the  introduction  of  the  insulating  microcracks  within 
the  framework  of  NIA. 

For  a  hexagonal  piezoelectric  solid  with  its  poling 
direction  along  thex3-axis,  it  fol lows  from  E q.  30  that 
the  effective  electroelastic  properties  of  the  micro- 
cracked,  piezoelectric  solid  pertaining  to  the  anti  plane 
deformation  and  in-plane  electric  fields  can  be  explic¬ 
itly  given  by 


r*  p*  c*  1 

c55  _  £15  _  fn  _ t _ 

C44  £15  £n  1  _|_  ^[200— ^sin(2flo)]  ’ 

^44  _  e24  _  e22  _  1 

C44  £15  £n  ]_  _|_  n *7[2flo+^sin(2fl())]  ’ 

45  —  e14  —  e25  —  e12  —  °- 


(32) 


For  a  hexagonal  piezoelectric  solid  with  its  poling 
direction  along  the  jc2-axis,  it  follows  from  Eq.  30  that 
the  effective  electroelastic  properties  of  the  micro- 
cracked,  piezoelectric  solid  pertaining  to  the  in-plane 
deformation  and  in-plane  electric  fields  can  be  explic¬ 
itly  given  by 


~C*U  Cf2  621  " 

« - 1 

1 

C11  C 12  621 

C12  C22  e22 

= 

C12  C22  e22 

_e21  e22  ~e22_ 

_  621  622  -622  _ 
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7tr]  [2$o  -  sin(20o)] 
40O 


7i r)  [20o  +  sin(20o)] 
+  40q 


r 0  °i 

0  0  0 

_  0  0  0_ 

ro  o  oi 


oi  i 

Ct  e 

0  1 

e  €  — 1 


'Ct 6 

e16 

-1 

C66 

eio 

.  e16 

-*uJ 

.  <?16 

-eii  _ 

7r ?7  [20o  -  sin(20o)] 
+  40o 

rr/7  [20o  +  sin(20o)] 
+  40q 


C*  /^*  n 

16  =  C26  ~  ell  -  e12  =  e26  -  e12  =  °> 


(33) 


where 
Re{Y}  = 


-  jl_ 

Cl 

0 

0 


0 


1 


e 


(34) 


has  been  defined  in  Eq.C12  by  Suo  et  al.  (1992).  The 
elements  in  Eq.34  have  to  be  determined  numerically 
for  the  specific  piezoelectric  material.  We  illustrate  in 
Figs.  3, 4, 5  the  variations  of  the  effective  electroelastic 
moduli  as  functions  of  0o  with  ri  =  0.2  for  a  micro- 
cracked  piezoelectric  material  B  aT i O 3  w i th  its  material 
properties  given  by  (Pan  2001) 


C 11  =  C33  =  166  x  109  N /m2, 

C13  =  77  x  109l\l/m2, 

C12  =  C23  =  78  x  109  N/m2, 

C22  =  162  x  109  N/m2, 

C44  =  C66  =  43  x  109  N/m2, 

C55  =  44.5  x  109  N /m2, 
e2i  =  e23  =  — 4.4C/m2, 
e?i  =  18.6C/m2, 

634  =  ei6  =  11.6  C/m2, 

en  =  e33  =  11.2  x  10_9C2/(Nm2), 

622  =  12.6  x  10_9C2/(Nm2). 


It  can  be  observed  from  the  three  figures  that  the 
degree  of  the  crack  alignment  characterized  by  Oo  has 
a  significant  influence  on  most  of  the  effective  elec¬ 
troelastic  moduli  except  C|6  due  to  the  fact  that  the 
ratio  Ct /Cl  =  1.0736  for  B aT i O 3  is  very  close  to 
unity.  We  also  observe  that  all  the  effective  electroelas¬ 
tic  moduli  are  monotonically  increasing  or  decreasing 


x1010(N/m2) 


Fig.  3  The  variations  of  the  four  effective  elastic  constants 
C{  1,  C22,  C{ 2  and  C6*6  as  functions  of  Oo  with  rj  =  0.2  for  a 
piezoelectric  material  BaTiO 3  containing  two-dimensional  insu¬ 
lating  mi  crocracks.  T  he  pol  i  ng  di  recti  on  of  the  pi  ezoel  ectri  c  sol  i  d 
is  along  the  jc2-axis 


(C/m2) 


Fig.  4  The  variations  of  the  three  effective  piezoelectric  con¬ 
stants  e\2,  ejj  and  e\b  as  functions  of  Oo  with  r/  =  0.2  for  a 
piezoelectric  material  BaTiO 3  containing  two-dimensional  insu- 
I  ati  ng  mi  crocracks.  T  he  pol  i  ng  di  recti  on  of  the  pi  ezoel  ectri  c  sol  i  d 
is  along  the *2-axis 


functions  of  Oo  at  a  fixed  crack  density  rj.  It  can  be 
clearly  seen  from  Fig.  3  that  the  effective  elastic  anisot¬ 
ropy  decreases  as  Oo  increases  from  zero  for  perfectly 
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Fig.  5  The  variations  of  the  two  effective  dielectric  constants 
e21  and  <?22  as  functions  of  9q  with  rj  =  0.2  for  a  piezoelectric 
material  BaTiCb  containing  two-dimensional  insulating  micro¬ 
cracks.  The  poling  direction  of  the  piezoelectric  solid  is  along 
thex2-axis 

aligned  microcracks  to  90°  for  randomly  oriented  mi¬ 
crocracks.  When  0q  =  90°,  Cj^  «  C|2  and  C|6  « 
\  (C^  -  C{2).  We  observe  from  Fig.  5  that  eft  >  en 
when  0o  <  40°.  This  observation  indicates  that  the 
piezoelectric  effect  plays  an  important  role  in  the  effec¬ 
tive  dielectric  properties  of  a  microcracked  piezoelec¬ 
tric  solid.  On  the  other  hand  our  results  show  that  some 
of  the  effective  piezoelectric  and  dielectric  moduli  at  a 
certain  fixed  0q  can  be  rather  complex  functions  of  rj. 
This  complex  phenomenon  is  demonstrated  in  Figs. 6 
and  7  for  e21  and  eft  as  functions  of  rj  for  four  differ¬ 
ent  values  of  0o  =  0,  n/8,  jt/4,  tt/2.  We  can  see  from 
Fig.  6  that  when  g'o  is  below  a  certain  value  (the  detailed 
results  show  that  this  value  is  0o  =  tt/3.03),  e21  can 
attain  its  minimum  value  at  an  extremely  low  crack 
density  rj  (for  example  when  0o  =  jt/8,  eft  attains  its 
minimum  value  of  -4.5851  C/m2  at  rj  =  0.04);  on  the 
other  hand  when  jt/3.03  <  0q  <  n/2,  eft  is  a  mono- 
tonically  increasing  function  of  rj.  Our  detailed  results 
also  show  that  depending  on  the  value  of  %  there  exist 
three  possible  distributions  of  eft:  (i)  a  monotonically 
increasing  function  of  rj  when  0o  =  0;  (ii)  a  function 
with  its  maximum  value  at  a  certain  value  of  ri  when 
0  <0o  <  zr/4.09  (for  example  as  shown  in  Fig.  7  when 
0o  =  7t/ 8,  eft  attains  its  maximum  value  of  11.86  x 
10_9C2/(l\lnr)  at  rj  =  0.436);  (iii)  a  monotonically 
decreasing  function  of  rj  when  0o  >  jt/4. 09.  We  add 


Fig.  6  The  variations  of  eft  as  a  function  of  o  at  four  differ¬ 
ent  values  of  e o  =  0,  n/8,  jt/4,  jt/2  for  a  piezoelectric  material 
BaTiCh  containing  two-dimensional  insulating  microcracks.  The 
poling  direction  of  the  piezoelectric  solid  is  along  thex2-axis 


Fig.  7  The  variations  of  eft  as  a  function  of  >i  at  four  differ¬ 
ent  values  of  $o  =  0,  jt/8,  jt/4,  jt/2  for  a  piezoelectric  material 
BaTi03  containing  two-dimensional  insulating  microcracks.  The 
poling  direction  of  the  piezoelectric  solid  is  along  thex2-axis 

that  this  unexpected  phenomenon  is  caused  by  the  cou¬ 
pling  effect  between  the  mechanical  field  and  the  elec¬ 
tric  field  (or  the  piezoelectric  effect).  Forexamplewhen 
there  exists  no  piezoelectric  effect,  eft  =  en  for  per¬ 
fectly  aligned  microcracks  (0o  =  0). 

3.2  A  piezoelectric  solid  containing  permeable 
microcracks 

The  effective  electroelastic  properties  of  a  piezoelec¬ 
tric  solid  containing  permeable  microcracks  (i.e.,  <p+  = 
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(p~,D+  =  D~  across  the  crack  surface,  see  Suo  et  al. 
1992)  can  be  simply  derived  as 

(C0*)-1  =  (C0)"1 

nrj  [2f?o  -  sin(20o)]  lTj  _iT 

I  7 -  J3  A-'ll  J3 


40O 

^•r?[20o  +  sin(20o)]T7'T  _iT 
+  40o  14  llJ4’ 


(35) 


where  the  3  x  3  symmetric  matrix  Ln  has  been  defined 
in  Eq.29  and 


Jb 


"1  0  0  0  0  0  0 

0  0  0  0  1  0  0 

0  0  0  1  0  0  0 


J4 


"0  0  0  0  1  0  0 

0  1  0  0  0  0  0 

0  0  1  0  0  0  0 


(36) 


We  see  from  Eq.35  that  the  thirty-two  components 
(12)  =  (21),  (13)  =  (31),  (24)  =  (42),  (34)  =  (43),  (ij)  = 
(ji),  (i  =  1-7,  j  =  6,7)  of  (C0*)-1  do  not  change  due  to 
the  introduction  of  the  permeable  microcracks  within 
the  framework  of  NIA.  For  a  hexagonal  piezoelectric 
solid  with  its  poling  direction  along  the*3-axis,  it  fol¬ 
lows  from  Eq.  35  that  the  effective  electroelastic  prop¬ 
erties  of  the  microcracked,  piezoelectric  solid  pertain¬ 
ing  to  the  anti  plane  deformation  and  in-plane  electric 
fields  can  be  explicitly  given  by 

C55  _  eh  _  _ l _ 

C44  ei5  1  _L  tt>?[26>0— sin(26>o)]  ’ 

x  -r  40() 

£11  _ 1  nk\i]  [20p  -  sin(20p)]  ^  ^ 

en  40o  +  nrj [20o  —  sin(20o)]  ~~ 

c*  p*  1 

c44  _  £24  _  _ t _ 

C44  ei5  ]_  _|_  jr,?[2c?o+si  n (2c?p)]  ’ 

£22  _  j  ,  zr/ci?7  [20q  +  sin(20o)]  >  ^ 

en  40o  +  nr)  [20o  +  sin(20o)] 

C45  =  =  el2  =  (37) 

where^i  =  e\b/(euC^)  <  1  is  the  electromechanical 
coupling  factor. 

We  observe  from  the  above  expression  that  only  the 
expressions  for  and  ef2  are  different  from  the  cor¬ 
responding  ones  in  Eq. 32  for  insulating  microcracks. 
When  0o  is  fixed,  and  e22  for  permeable  micro¬ 
cracks  are  increasing  functions  of  +  whereas  those  for 
insulating  microcracks  are  decreasing  functions  of  rj. 

Similarly,  for  a  hexagonal  piezoelectric  solid  with 
its  poling  direction  along  the  x2-axis,  it  follows  from 
Eq.  35  that  the  effective  electroelastic  properties  of  the 


microcracked  piezoelectric  solid  pertaining  to  the  in¬ 
plane  deformation  and  in-plane  electric  fields  can  be 
explicitly  given  by 
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(38) 


where  ki  =  e2/(eCT). 


3.3  A  piezoelectric  solid  containing  conducting 
microcracks 


The  eff ecti ve el ectroel asti c  properti es  of  a  pi ezoel ectri c 
solid  containing  conducting  microcracks  (M  cM  eeking 
1987;  Suo  1993)  can  also  be  similarly  derived.  The 
energy  balance  relationship  for  two-dimensional  prob¬ 
lems  can  be  expressed  in  the  following  form  when  the 
piezoelectric  solid  containing  conducting  microcracks 
is  subjected  to  the  homogeneous  stresses  crP.  and  the 
homogeneous  electric  fields  Ef  (Suo  1993) 


1  — 1 X- 0  ^  /  x.  0.7,  V,  0  \  — 1 X- 0 
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(39) 


where 
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[u]  =  [[«!]  [M2]  [MB]  tm3 


(41) 
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(44) 


In  E q. 41,  £  is  defined  such  that  D\  =  q,  £>2  = 
-£1.  In  Eqs.42  and  43,  q,  Ai;-  and  frj  are  the  elas¬ 
tic,  piezoelectric  and  dielectric  constants  of  the  undam¬ 


aged  solid,  while  C*j 


h*j  and  pfj  arethecorresponding 


(unknown)  constants  for  the  damaged  solid.  Here  we 
have  adopted  the  material  constant  notations  in  Suo 
(1993).  When  ignoring  crack  interactions,  the  jump  in 
displacement  and  £  across  the  crack  surfaces  can  be 
simply  given  by  (Suo  1993) 

[u]  =  2a/g2  -*2L_1t°,  (|x|  <  a)  (45) 

where  L  is  a  4  x  4  real  and  symmetric,  positive  definite 
matrix. 

Substituting  the  above  results  into  Eq.  22,  we  finally 
arrive  at  a  concise  expression  for  the  effective  elec¬ 
troelastic  properties  of  the  microcracked  piezoelectric 
solid 
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(47) 


It  is  apparent  that  the  energy  density  function  always 
increases  through  the  introduction  of  conducting  mi¬ 
crocracks,  and  that  C°*  is  positive  definite  in  view  of 
the  fact  that  both  C°  and  L  are  positive  definite.  I  nter- 
esti  ng  I  y,  th i  s  resu 1 1  i  s  very  si  m  i  I  ar  to  the  si  tu ati  0 n  w  h en 
the  matrix  is  purely  elastic.  We  stress  that  the  convex 
energy  density  function  used  here  for  conducting  cracks 
is  different  than  the  non-convex  electric  enthalpy  used 
for  insulating  cracks. 

For  a  hexagonal  piezoelectric  solid  with  its  poling 
direction  along  thex3-axis,  it  fol lows  from  Eq.  46  that 
the  effective  electroelastic  properties  of  the  micro- 
cracked  piezoelectric  solid  pertaining  to  the  anti  plane 
deformation  and  in-plane  electric  fields  can  be  explic¬ 
itly  given  by 
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where C44P11  >  hh.  A pparently  c|4,  C|5,  p{v  p22  > 

0  and  ^55/^33  >  ^15’  £'44  022  >  *24' 

For  a  hexagonal  piezoelectric  solid  with  its  poling 
direction  along  thex2-axis,  it  follows  from  Eq.  46  that 
the  effective  electroelastic  properties  of  the  micro- 
cracked,  piezoelectric  solid  pertaining  to  the  in-plane 
deformation  and  in-plane  electric  fields  can  be  explic¬ 
itly  given  by 
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nr]  [20o  -  sin(20o)] 
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nrj  [20o  +  sin(20o)] 
40O 


Si  d 
d  e 


C  f6  =  C2*6  =  h*u  =  h*u  =  h*26  =  =  0,  (49) 

where 


Re{Y)  = 


Si  0  d 
0  S30 
d  0  f 


(50) 


has  been  defined  by  Suo  (1993).  Similar  to  the 
case  of  insulating  cracks,  the  elements  in  E q. 50 
have  to  be  computed  numerically.  Apparently  both 


r*  r1*  u* 

C11  C12  n21 
r '*  n*  h* 
c12  c22  n22 

-h21  h22  &22- 


and 


'^66*16~ 
^i6  Phi 


are  positive  definite. 


which  remains  valid  when  p  «;  1.  We  see  from  the 
above  expression  that  the  structure  of  k  is  exactly  the 
same  as  that  of  L  described  in  Eq.29. 


3.4.1  Remark 


Fan  and  Sze  (2001)  derived  a  closed-form  expression  in 
(2.9)  forthespring  constant*  for  two-dimensional  non¬ 
interacting  cracks.  In  fact,  the  expression  for  k  in  three- 
dimensions  for  penny-shaped,  non-interacting  cracks 
can  also  be  similarly  derived  as 


k  = 


3nC  1 

8 a  p2  ’ 


(54) 


which  matches  quite  well  the  finite  element  results  in 
Table  2  in  Fan  and  Sze  (2001)  up  to  p  =  0.5. 


3.4  A  special  distribution  of  the  insulating 
microcracks 

Next  we  discuss  a  special  distribution  of  the  micro¬ 
cracks:  all  the  insulating  microcracks  with  common 
half-length  a  are  distributed  along  the  xi-axis  in  a 
homogeneous  piezoelectric  solid  subjected  to  the 
homogeneous  stresses  <tP  and  the  homogeneous 
electric  displacements  £>9.  As  pointed  out  by  Fan  and 
Sze  (2001),  we  may  merge  the  microcracks  along  the 
xi-axi s  i  nto  a  conti  nuously  damaged  i  nterface  such  that 

(t2)=k<[u]>,  on  x2  =  0  (51) 

where  (*)  stands  for  the  average  over  an  area  of  a  scale 
much  greater  than  the  dimension  of  the  microcracks, 
t2  =  [cti2  cr 22  <732  D2V  and  [u]  being  defined  by 
Eq.24.  In  addition  (t2>  =  t°  =  [ct°2  ct2°2  °32  DVY  ■ 
E  quati  on  51  can  be  termed  the  general  ized,  spri  ng-ty  pe, 
imperfect  interface  model  with  the  4  x  4  symmetric 
matrix  k  being  the  generalized  stiffness  matrix  to  be 
determined.  When  ignoring  the  crack  interactions  and 
noticing  Eqs.  12  and  28  with  0  =  0,  we  can  easily 
obtain  the  following 

<[u]>  =  (52) 

where  p  =  Ma/L  with  M  being  the  number  of  mi¬ 
crocracks  within  the  sample  of  half-length  L  taken  out 
of  the  jq-axis.  Comparison  of  Eq.52  with  E q. 51  will 
immediately  lead  to  the  foil  owing  closed-form  expres¬ 
sion  of  k  as 

k=— L,  (53) 

nap 


4  An  approximation  of  the  GSCM 


We  point  out  that  the  problem  in  Sect.  2  can  also  be 
discussed  within  theframework  of  the  general  ized  self- 
consi  stent  method(GSCM),whichapproxi  matel  y  takes 
into  account  crack  interaction  (Aboudi  and  Benveniste 
1987;  Santareetal.  1995).  Herewecan  adopt  a  simpli¬ 
fied  version  of  the  GSCM  for  a  microcracked,  elasti¬ 
cally  anisotropic  sol  id:  first  calculate  the  uniform  stress 
field  <7°  withinanintactdrcM/arinclusion  with  undam¬ 
aged  material  properties  surrounded  by  the  effective 
medium  with  (unknown)  material  properties  subjected 
to  remote  uniform  stresses  <t?.  (Ting  1996);  second 
solve  the  probl  em  of  a  crack  in  an  infinite  homogeneous 
material  with  undamaged  material  properties  subjected 
to  remote  uniform  loading  o-P  obtained  in  step  1.  After 
adopting  the  above,  simplified  version  of  the  GSCM  , 
the  crack  opening  displacement  can  be  finally 
obtained  as 

[u]  =  2 V a2  —  x2L-1t,  (|x|  <  a)  (55) 

where 

t  =  L(Eit°  +  E2t§)  COS  0  -  L(FitJ  +  F2t§)  Sin  0, 

(56) 
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In  the  above  expressions  the  quantities  attached  with 
the  superscript  *  are  those  pertaining  to  the  effective 
medium.  H  ere  we  point  out  that  the  above  derivations 
are  based  on  the  el  I  i  pti  c  i  ncl  usi  on  and  crack  sol  uti  ons  i  n 
Ting  (1996)  and  more  detailed  structures  and  identities 
of  Ni,  N2,  N3,S,HandL,  all  of  which  can  be  explicitly 
expressed  in  terms  of  the  reduced  elastic  compliances, 
can  be  found  in  Chapter  6  of  Ting  (1996).  The  Bar¬ 
nett- Lothe  tensor,  S,  should  not  be  confused  with  the 
elastic  compliance.  It  can  be  easily  verified  from  the 
above  expressions  that  t  =  cos  9  -  tj  sin  0  when 
the  circular  inclusion  and  the  surrounding  matrix  pos¬ 
sess  exactly  the  same  material  properties.  Substituting 
the  above  results  into  the  energy  balance  expression, 
Eq.3,  we  finally  arrive  at  the  following  set  of  fifteen 
coupled,  nonlinear  equations  for  the  fifteen  unknown 
effective  reduced  elastic  compliances: 
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(63) 

where  E,V,  £ (2) ,  F;[P  and  F®,  all  of  which  arefunc- 

i  J  l  J  i  J  i  J 

tions  of  the  unknown  effective  reduced  elastic  com¬ 
pliances,  are  the  components  of  four  3x3  matrices 
Ei,  E2,  Fi  and  F2  respectively. 

Equations61-63  can  be  solved  through  iteration 
(Santare  et  al.  1995;  Wang  et  al.  2009).  During  spe¬ 
cific  iteration  for  the  unknowns,  it  is  more  convenient 
to  adopt  the  following  equivalent  matrix  expression 
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We  observe  that  when  using  the  GSCM ,  the  fifteen 
effective,  reduced  elastic  compliances  will,  in  general, 
be  different  from  the  corresponding  compliances  for 
the  undamaged  material.  Thus  the  result  that  Sn  = 

S'i2’  ^14  =  ^25  =  ^25’  ^45  =  “^45  Eg.  14)  is 

only  valid  for  NIA. 


5  Effective  thermal  conductivities  of  a  solid 
containing  insulated  penny-shaped  cracks 


I  n  the  above,  the  crack  orientation  distribution  function 
4>(e)  has  been  introduced  for  two-dimensional  cracks. 
In  fact,  the  crack  orientation  distribution  function  can 
also  be  adopted  to  describe  three-dimensional  cracks. 
As  an  illustration,  we  consider  the  effective  thermal 
conductivities  of  a  solid  containing  insulated  penny¬ 
shaped  cracks  with  an  arbitrary  degree  of  alignment. 
Within  the  framework  of  steady-state  heat  conduction, 
we  consider  insulated  penny-shaped  cracks  of  radius  a 
with  an  arbitrary  degree  of  alignment,  distributed  in  a 
homogeneous  and  isotropic  medium  with  thermal  con¬ 
ductivity  k.  The  orientation  of  the  penny-shaped  crack 
can  be  described  by  the  azimuth  angle  9  (\9\  <  n/2) 
and  the  zenith  angle  </>  (0  <  <p  <  jt)  which  define 
the  direction  of  the  unit  normal  of  the  crack  surface  as 
shown  in  Fig.  8  (Yang  and  Turner  2003).  The  crack  ori¬ 
entation  distribution  function  (9,  <fr)  on  a  unit  sphere 
can  be  given  by 


<P(9,<p ) 

AQ  ■ — .  |0|  <  and  j  —  <po  <  <p  <  j  +  <po 
=  •  4$o  si  n  cpo 

0,  else 

(66) 

where  9q,  <po  <  n/2.  By  adopting  the  energy  method 
and  ignoring  the  crack  interactions,  closed-form 
expressions  for  the  effective  thermal  conductivities  can 
be  derived  as 
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*3*3  =  - z-T!— ,  (67) 

i  ,  8// sin^  ifio 

1  H - 5 - 

where  rj  =  Ma3/  V  is  the  crack  density  parameter  and 
k*j  =  0,  (/  ^  /),  which  implies  that  the  three  prin¬ 
cipal  directions  of  the  effective  medium  are  along  the 
xi,  X2  and  X3  axes.  The  above  derivation  uses  the  fol¬ 
lowing  expression  for  the  temperature  jump  [r]  across 
the  insulated  crack  surfaces  (Barber  1975) 

4qo  —  r2 


[T]  =  — 


7TK 


(0  <  r  <  a) 


(68) 


where  qo  is  the  uniform  heat  flow  prescribed  on  the 
crack  face. 

If  all  the  microcracks  are  randomly  oriented,  i.e., 
(p0  =  e0  =  tt/2,  then  we  have 

Kii  =  *22  =  ^33  =  rrr 


1  +  5)7 


(69) 


which  means  that  the  effective  property  is  still  isotro¬ 
pic. 

If  the  unit  normals  of  all  the  microcracks  lie  in  the 

xi  -  X2  plane,  i.e.,  y>o  =  0,  then  we  have 

*  K 
ku  = 


i  ,  2r,[20o+sm(2eo)]  ’ 
i+  355 


K22  = 


i  ,  2»,[20o-sin(20o)]  ’ 
i+  Wo 


*3*3  =  *.  (70) 

which  indicates  that  the  effective  thermal  conductiv¬ 
ity  /cf3  is  not  influenced  by  the  existence  of  the  mi¬ 
crocracks.  Furthermore,  if  the  unit  normals,  which  lie 
in  the  xi  -  X2  plane,  have  a  random  distribution,  i.e., 
9o  =  it/ 2  (or  the  so  called  uniaxially  aligned  cracks, 
seeYang  and  Turner  2003),  then  itfollows  from  Eq.  70 
that 

=  k22  =  Kh  =  (71> 

1+3  ^ 
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which  indicates  that  the  effective  property  is  trans¬ 
versely  isotropic  with  the  x3-axis  taken  as  the  uniaxial 
symmetry  axis.  Alternatively,  if  the  unit  normals  lying 
in  the  xi  -  X2  plane  are  along  thexi-axis,  i.e.,  6>o  =  0 
(i.e,  the  so-called  perfectly  aligned  cracks,  see  Yang 
and  Turner  2005),  then  it  follows  from  Eq.  70  that 

*n  =  i  ~ 8  ’  Kh  =  *33  =  *>  <72> 

1  +  §/? 

which  indicates  that  the  effective  property  is  also  trans¬ 
versely  isotropic,  but  now  with  the  jq-axis  as  the 
uniaxial  symmetry  axis. 


6  Conclusions 

We  first  derived  in  Eq.14  the  closed-form  expressions 
for  the  fifteen  effective  reduced  elastic  compliances  of 
a  microcracked,  elastically  anisotropic  solid  by  using 
NIA.  Then  we  further  obtained  in  Eqs.30,  35  and  46 
the  concise  expressions  for  the  effective  electroelastic 
properties  of  a  piezoelectric  solid  containing  insulat¬ 
ing,  permeable  or  conducting  microcracks  by  extend¬ 
ing  the  N I A  for  purel y  el  asti  c  materi  al  s.  W e  al  so  derived 
a  set  of  coupled  nonl  inear  Eq.  64  for  the  unknown  effec¬ 
tive  reduced  elastic  compliances  by  using  the  GSCM  . 
Finally  the  crack  orientation  distribution  function  was 
introduced  to  investigate  the  effective  thermal  conduc¬ 
tivities  of  a  solid  containing  insulated  penny-shaped 
cracks  with  an  arbitrary  degree  of  alignment,  and  the 
effective  thermal  conductivities  were  derived  in  Eq.67. 
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